source: trunk/GSASIIindex.py @ 1590

Last change on this file since 1590 was 1587, checked in by vondreele, 11 years ago

use successive approx. to get d from TOF
allow sorting on a,b,c,alp,bet,gam, vol & m20 in cell indexed table
put in derivatives for fitCell inside indexing routine - better accuracy & speed
add a couple more orthos to SS table
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 36.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII cell indexing program: variation on that of A. Coehlo
3#   includes cell refinement from peak positions (not zero as yet)
4########### SVN repository information ###################
5# $Date: 2014-11-27 15:34:45 +0000 (Thu, 27 Nov 2014) $
6# $Author: toby $
7# $Revision: 1587 $
8# $URL: trunk/GSASIIindex.py $
9# $Id: GSASIIindex.py 1587 2014-11-27 15:34:45Z toby $
10########### SVN repository information ###################
11'''
12*GSASIIindex: Cell Indexing Module*
13===================================
14
15Cell indexing program: variation on that of A. Coehlo
16includes cell refinement from peak positions (not zero as yet)
17
18This needs a bit of refactoring to remove the little bit of GUI
19code referencing wx
20'''
21
22import math
23import wx
24import time
25import numpy as np
26import numpy.linalg as nl
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 1587 $")
29import GSASIIlattice as G2lat
30import GSASIIpwd as G2pwd
31import scipy.optimize as so
32
33# trig functions in degrees
34sind = lambda x: math.sin(x*math.pi/180.)
35asind = lambda x: 180.*math.asin(x)/math.pi
36tand = lambda x: math.tan(x*math.pi/180.)
37atand = lambda x: 180.*math.atan(x)/math.pi
38atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
39cosd = lambda x: math.cos(x*math.pi/180.)
40acosd = lambda x: 180.*math.acos(x)/math.pi
41rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
42#numpy versions
43npsind = lambda x: np.sin(x*np.pi/180.)
44npasind = lambda x: 180.*np.arcsin(x)/math.pi
45npcosd = lambda x: np.cos(x*math.pi/180.)
46nptand = lambda x: np.tan(x*math.pi/180.)
47npatand = lambda x: 180.*np.arctan(x)/np.pi
48npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
49rpd = np.pi/180.
50   
51def scaleAbyV(A,V):
52    'needs a doc string'
53    v = G2lat.calc_V(A)
54    scale = math.exp(math.log(v/V)/3.)**2
55    for i in range(6):
56        A[i] *= scale
57   
58def ranaxis(dmin,dmax):
59    'needs a doc string'
60    import random as rand
61    return rand.random()*(dmax-dmin)+dmin
62   
63def ran2axis(k,N):
64    'needs a doc string'
65    import random as rand
66    T = 1.5+0.49*k/N
67#    B = 0.99-0.49*k/N
68#    B = 0.99-0.049*k/N
69    B = 0.99-0.149*k/N
70    R = (T-B)*rand.random()+B
71    return R
72   
73#def ranNaxis(k,N):
74#    import random as rand
75#    T = 1.0+1.0*k/N
76#    B = 1.0-1.0*k/N
77#    R = (T-B)*rand.random()+B
78#    return R
79   
80def ranAbyV(Bravais,dmin,dmax,V):
81    'needs a doc string'
82    cell = [0,0,0,0,0,0]
83    bad = True
84    while bad:
85        bad = False
86        cell = rancell(Bravais,dmin,dmax)
87        G,g = G2lat.cell2Gmat(cell)
88        A = G2lat.Gmat2A(G)
89        if G2lat.calc_rVsq(A) < 1:
90            scaleAbyV(A,V)
91            cell = G2lat.A2cell(A)
92            for i in range(3):
93                bad |= cell[i] < dmin
94    return A
95   
96def ranAbyR(Bravais,A,k,N,ranFunc):
97    'needs a doc string'
98    R = ranFunc(k,N)
99    if Bravais in [0,1,2]:          #cubic - not used
100        A[0] = A[1] = A[2] = A[0]*R
101        A[3] = A[4] = A[5] = 0.
102    elif Bravais in [3,4]:          #hexagonal/trigonal
103        A[0] = A[1] = A[3] = A[0]*R
104        A[2] *= R
105        A[4] = A[5] = 0.       
106    elif Bravais in [5,6]:          #tetragonal
107        A[0] = A[1] = A[0]*R
108        A[2] *= R
109        A[3] = A[4] = A[5] = 0.       
110    elif Bravais in [7,8,9,10]:     #orthorhombic
111        A[0] *= R
112        A[1] *= R
113        A[2] *= R
114        A[3] = A[4] = A[5] = 0.       
115    elif Bravais in [11,12]:        #monoclinic
116        A[0] *= R
117        A[1] *= R
118        A[2] *= R
119        A[4] *= R
120        A[3] = A[5] = 0.       
121    else:                           #triclinic
122        A[0] *= R
123        A[1] *= R
124        A[2] *= R
125        A[3] *= R
126        A[4] *= R
127        A[5] *= R
128    return A
129   
130def rancell(Bravais,dmin,dmax):
131    'needs a doc string'
132    if Bravais in [0,1,2]:          #cubic
133        a = b = c = ranaxis(dmin,dmax)
134        alp = bet = gam = 90
135    elif Bravais in [3,4]:          #hexagonal/trigonal
136        a = b = ranaxis(dmin,dmax)
137        c = ranaxis(dmin,dmax)
138        alp = bet =  90
139        gam = 120
140    elif Bravais in [5,6]:          #tetragonal
141        a = b = ranaxis(dmin,dmax)
142        c = ranaxis(dmin,dmax)
143        alp = bet = gam = 90
144    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
145        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
146        abc.sort()
147        a = abc[0]
148        b = abc[1]
149        c = abc[2]
150        alp = bet = gam = 90
151    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
152        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
153        ac.sort()
154        a = ac[0]
155        b = ranaxis(dmin,dmax)
156        c = ac[1]
157        alp = gam = 90
158        bet = ranaxis(90.,130.)
159    else:                           #triclinic - a<b<c convention
160        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
161        abc.sort()
162        a = abc[0]
163        b = abc[1]
164        c = abc[2]
165        r = 0.5*b/c
166        alp = ranaxis(acosd(r),acosd(-r))
167        r = 0.5*a/c
168        bet = ranaxis(acosd(r),acosd(-r))
169        r = 0.5*a/b
170        gam = ranaxis(acosd(r),acosd(-r)) 
171    return [a,b,c,alp,bet,gam]
172   
173def calc_M20(peaks,HKL):
174    'needs a doc string'
175    diff = 0
176    X20 = 0
177    for Nobs20,peak in enumerate(peaks):
178        if peak[3]:
179            Qobs = 1.0/peak[7]**2
180            Qcalc = 1.0/peak[8]**2
181            diff += abs(Qobs-Qcalc)
182        elif peak[2]:
183            X20 += 1
184        if Nobs20 == 19: 
185            d20 = peak[7]
186            break
187    else:
188        d20 = peak[7]
189        Nobs20 = len(peaks)
190    for N20,hkl in enumerate(HKL):
191        if hkl[3] < d20:
192            break               
193    eta = diff/Nobs20
194    Q20 = 1.0/d20**2
195    if diff:
196        M20 = Q20/(2.0*diff)
197    else:
198        M20 = 0
199    M20 /= (1.+X20)
200    return M20,X20
201   
202def calc_M20SS(peaks,HKL):
203    'needs a doc string'
204    diff = 0
205    X20 = 0
206    for Nobs20,peak in enumerate(peaks):
207        if peak[3]:
208            Qobs = 1.0/peak[8]**2
209            Qcalc = 1.0/peak[9]**2
210            diff += abs(Qobs-Qcalc)
211        elif peak[2]:
212            X20 += 1
213        if Nobs20 == 19: 
214            d20 = peak[8]
215            break
216    else:
217        d20 = peak[8]
218        Nobs20 = len(peaks)
219    for N20,hkl in enumerate(HKL):
220        if hkl[4] < d20:
221            break               
222    eta = diff/Nobs20
223    Q20 = 1.0/d20**2
224    if diff:
225        M20 = Q20/(2.0*diff)
226    else:
227        M20 = 0
228    M20 /= (1.+X20)
229    return M20,X20
230   
231def sortM20(cells):
232    'needs a doc string'
233    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
234    #sort highest M20 1st
235    T = []
236    for i,M in enumerate(cells):
237        T.append((M[0],i))
238    D = dict(zip(T,cells))
239    T.sort()
240    T.reverse()
241    X = []
242    for key in T:
243        X.append(D[key])
244    return X
245               
246def sortCells(cells,col):
247    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam,volume
248    #sort smallest a,b,c,alpha,beta,gamma or volume 1st
249    T = []
250    for i,M in enumerate(cells):
251        T.append((M[col],i))
252    D = dict(zip(T,cells))
253    T.sort()
254    X = []
255    for key in T:
256        X.append(D[key])
257    return X
258   
259def findMV(peaks,HKL,ssopt):
260#    import basinhopping as bh
261    print ssopt
262    return ssopt['ModVec']
263               
264def IndexPeaks(peaks,HKL):
265    'needs a doc string'
266    import bisect
267    N = len(HKL)
268    if N == 0: return False,peaks
269    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
270    hklds.sort()                                        # ascending sort - upper bound at end
271    hklmax = [0,0,0]
272    for ipk,peak in enumerate(peaks):
273        if peak[2]:
274            peak[4:7] = [0,0,0]                           #clear old indexing
275            peak[8] = 0.
276            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
277            dm = peak[-2]-hklds[i-1]                         # peak to neighbor hkls in list
278            dp = hklds[i]-peak[-2]
279            pos = N-i                                       # reverse the order
280            if dp > dm: pos += 1                            # closer to upper than lower
281            if pos >= N:
282                print pos,N
283                break
284            hkl = HKL[pos]                                 # put in hkl
285            if hkl[-1] >= 0:                                 # peak already assigned - test if this one better
286                opeak = peaks[hkl[-1]]
287                dold = abs(opeak[-2]-hkl[3])
288                dnew = min(dm,dp)
289                if dold > dnew:                             # new better - zero out old
290                    opeak[4:7] = [0,0,0]
291                    opeak[8] = 0.
292                else:                                       # old better - do nothing
293                    continue               
294            hkl[-1] = ipk
295            peak[4:7] = hkl[:3]
296            peak[8] = hkl[3]                                # fill in d-calc
297    for peak in peaks:
298        peak[3] = False
299        if peak[2]:
300            if peak[-1] > 0.:
301                for j in range(3):
302                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
303                peak[3] = True
304    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
305        return True,peaks
306    else:
307        return False,peaks  #nothing indexed!
308       
309def IndexSSPeaks(peaks,HKL):
310    'needs a doc string'
311    import bisect
312    N = len(HKL)
313    if N == 0: return False,peaks
314    if len(peaks[0]) == 9:      #add m column if missing
315        for peak in peaks:
316            peak.insert(7,0)
317    hklds = list(np.array(HKL).T[4])+[1000.0,0.0,]
318    hklds.sort()                                        # ascending sort - upper bound at end
319    hklmax = [0,0,0,0]
320    for ipk,peak in enumerate(peaks):
321        if peak[2]: #Use
322            peak[4:8] = [0,0,0,0]                           #clear old indexing
323            peak[9] = 0.
324            i = bisect.bisect_right(hklds,peak[8])          # find peak position in hkl list
325            dm = peak[8]-hklds[i-1]                         # peak to neighbor hkls in list
326            dp = hklds[i]-peak[8]
327            pos = N-i                                       # reverse the order
328            if dp > dm: pos += 1                            # closer to upper than lower
329            if pos >= N:
330                print pos,N
331                break
332            hkl = HKL[pos]                                 # put in hkl
333            if hkl[-1] >= 0:                                 # peak already assigned - test if this one better
334                opeak = peaks[hkl[-1]]
335                dold = abs(opeak[-2]-hkl[4])
336                dnew = min(dm,dp)
337                if dold > dnew:                             # new better - zero out old
338                    opeak[4:8] = [0,0,0,0]
339                    opeak[9] = 0.
340                else:                                       # old better - do nothing
341                    continue
342            hkl[-1] = ipk
343            peak[4:8] = hkl[:4]
344            peak[9] = hkl[4]                                # fill in d-calc
345    for peak in peaks:
346        peak[3] = False
347        if peak[2]:
348            if peak[-1] > 0.:
349                for j in range(4):
350                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
351                peak[3] = True
352    if hklmax[0]*hklmax[1]*hklmax[2]*hklmax[3] > 0:
353        return True,peaks
354    else:
355        return False,peaks  #nothing indexed!
356       
357def Values2A(ibrav,values):
358    'needs a doc string'
359    if ibrav in [0,1,2]:
360        return [values[0],values[0],values[0],0,0,0]
361    elif ibrav in [3,4]:
362        return [values[0],values[0],values[1],values[0],0,0]
363    elif ibrav in [5,6]:
364        return [values[0],values[0],values[1],0,0,0]
365    elif ibrav in [7,8,9,10]:
366        return [values[0],values[1],values[2],0,0,0]
367    elif ibrav in [11,12]:
368        return [values[0],values[1],values[2],0,values[3],0]
369    else:
370        return list(values[:6])
371       
372def A2values(ibrav,A):
373    'needs a doc string'
374    if ibrav in [0,1,2]:
375        return [A[0],]
376    elif ibrav in [3,4,5,6]:
377        return [A[0],A[2]]
378    elif ibrav in [7,8,9,10]:
379        return [A[0],A[1],A[2]]
380    elif ibrav in [11,12]:
381        return [A[0],A[1],A[2],A[4]]
382    else:
383        return A
384       
385def Values2Vec(ibrav,vec,Vref,val):
386    if ibrav in [3,4,5,6]:
387        Nskip = 2
388    elif ibrav in [7,8,9,10]:
389        Nskip = 3
390    elif ibrav in [11,12]:
391        Nskip = 4
392    else:
393        Nskip = 6
394    Vec = []
395    i = 0
396    for j,r in enumerate(Vref):
397        if r:
398            Vec.append(val[i+Nskip])
399            i += 1
400        else:
401            Vec.append(vec[j])
402    return np.array(Vec) 
403
404def FitHKL(ibrav,peaks,A,Pwr):
405    'needs a doc string'
406               
407    def errFit(values,ibrav,d,H,Pwr):
408        A = Values2A(ibrav,values)
409        Qo = 1./d**2
410        Qc = G2lat.calc_rDsq(H,A)
411        return (Qo-Qc)*d**Pwr
412       
413    def dervFit(values,ibrav,d,H,Pwr):
414        if ibrav in [0,1,2]:
415            derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
416        elif ibrav in [3,4,]:
417            derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
418        elif ibrav in [5,6]:
419            derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
420        elif ibrav in [7,8,9,10]:
421            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
422        elif ibrav in [11,12]:
423            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
424        else:
425            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
426        derv = -np.array(derv)
427        return (derv*d**Pwr).T
428   
429    Peaks = np.array(peaks).T
430    values = A2values(ibrav,A)
431    result = so.leastsq(errFit,values,Dfun=dervFit,full_output=True,ftol=0.000001,
432        args=(ibrav,Peaks[7],Peaks[4:7],Pwr))
433    A = Values2A(ibrav,result[0])
434    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
435           
436def errFitZ(values,ibrav,d,H,tth,wave,Z,Zref):
437    Zero = Z
438    if Zref:   
439        Zero = values[-1]
440    A = Values2A(ibrav,values)
441    Qo = 1./d**2
442    Qc = G2lat.calc_rDsqZ(H,A,Zero,tth,wave)
443    return (Qo-Qc)
444   
445def dervFitZ(values,ibrav,d,H,tth,wave,Z,Zref):
446    if ibrav in [0,1,2]:
447        derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
448    elif ibrav in [3,4,]:
449        derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
450    elif ibrav in [5,6]:
451        derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
452    elif ibrav in [7,8,9,10]:
453        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
454    elif ibrav in [11,12]:
455        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
456    else:
457        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
458    if Zref:
459        derv.append(npsind(tth)*2.0*rpd/wave**2)
460    derv = -np.array(derv)
461    return derv.T
462   
463def FitHKLZ(wave,ibrav,peaks,A,Z,Zref):
464    'needs a doc string'
465   
466    Peaks = np.array(peaks).T   
467    values = A2values(ibrav,A)
468    if Zref:
469        values.append(Z)
470    result = so.leastsq(errFitZ,values,Dfun=dervFitZ,full_output=True,ftol=0.0001,
471        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref))
472    A = Values2A(ibrav,result[0][:6])
473    if Zref:
474        Z = result[0][-1]
475    chisq = np.sum(errFitZ(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref)**2)
476    return True,chisq,A,Z,result
477   
478def errFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
479    Zero = Z
480    if Zref:   
481        Zero = values[-1]
482    A = Values2A(ibrav,values)
483    Vec = Values2Vec(ibrav,vec,Vref,values)
484    Qo = 1./d**2
485    Qc = G2lat.calc_rDsqZSS(H,A,Vec,Zero,tth,wave)
486    return (Qo-Qc)
487   
488def dervFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
489    A = Values2A(ibrav,values)
490    Vec = Values2Vec(ibrav,vec,Vref,values)
491    HM = H[:3]+(H[3][:,np.newaxis]*Vec).T
492    if ibrav in [3,4,]:
493        derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]]
494    elif ibrav in [5,6]:
495        derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]]
496    elif ibrav in [7,8,9,10]:
497        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]]
498    elif ibrav in [11,12]:
499        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]]
500    else:
501        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[1],HM[0]*HM[2],HM[1]*HM[2]]
502    if Vref[0]:
503        derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3])
504    if Vref[1]:
505        derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3])
506    if Vref[2]:
507        derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3])   
508    if Zref:
509        derv.append(npsind(tth)*2.0*rpd/wave**2)
510    derv = -np.array(derv)
511    return derv.T
512   
513def FitHKLZSS(wave,ibrav,peaks,A,V,Vref,Z,Zref):
514    'needs a doc string'
515   
516    Peaks = np.array(peaks).T   
517    values = A2values(ibrav,A)
518    for v,r in zip(V,Vref):
519        if r:
520            values.append(v)
521    if Zref:
522        values.append(Z)
523    result = so.leastsq(errFitZSS,values,Dfun=dervFitZSS,full_output=True,ftol=1.e-6,
524        args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,V,Vref,Z,Zref))
525    A = Values2A(ibrav,result[0])
526    Vec = Values2Vec(ibrav,V,Vref,result[0])
527    if Zref:
528        Z = result[0][-1]
529    chisq = np.sum(errFitZSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,Vec,Vref,Z,Zref)**2) 
530    return True,chisq,A,Vec,Z,result
531   
532def errFitT(values,ibrav,d,H,tof,difC,Z,Zref):
533    Zero = Z
534    if Zref:   
535        Zero = values[-1]
536    A = Values2A(ibrav,values)
537    Qo = 1./d**2
538    Qc = G2lat.calc_rDsqT(H,A,Zero,tof,difC)
539    return (Qo-Qc)
540   
541def dervFitT(values,ibrav,d,H,tof,difC,Z,Zref):
542    if ibrav in [0,1,2]:
543        derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
544    elif ibrav in [3,4,]:
545        derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
546    elif ibrav in [5,6]:
547        derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
548    elif ibrav in [7,8,9,10]:
549        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
550    elif ibrav in [11,12]:
551        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
552    else:
553        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
554    if Zref:
555        derv.append(np.ones_like(d)/difC)
556    derv = -np.array(derv)
557    return derv.T
558   
559def FitHKLT(difC,ibrav,peaks,A,Z,Zref):
560    'needs a doc string'
561   
562    Peaks = np.array(peaks).T   
563    values = A2values(ibrav,A)
564    if Zref:
565        values.append(Z)
566    result = so.leastsq(errFitT,values,Dfun=dervFitT,full_output=True,ftol=0.0001,
567        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref))
568    A = Values2A(ibrav,result[0])
569    if Zref:
570        Z = result[0][-1]
571    chisq = np.sum(errFitT(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref)**2)
572    return True,chisq,A,Z,result
573               
574def rotOrthoA(A):
575    'needs a doc string'
576    return [A[1],A[2],A[0],0,0,0]
577   
578def swapMonoA(A):
579    'needs a doc string'
580    return [A[2],A[1],A[0],0,A[4],0]
581   
582def oddPeak(indx,peaks):
583    'needs a doc string'
584    noOdd = True
585    for peak in peaks:
586        H = peak[4:7]
587        if H[indx] % 2:
588            noOdd = False
589    return noOdd
590   
591def halfCell(ibrav,A,peaks):
592    'needs a doc string'
593    if ibrav in [0,1,2]:
594        if oddPeak(0,peaks):
595            A[0] *= 2
596            A[1] = A[2] = A[0]
597    elif ibrav in [3,4,5,6]:
598        if oddPeak(0,peaks):
599            A[0] *= 2
600            A[1] = A[0]
601        if oddPeak(2,peaks):
602            A[2] *=2
603    else:
604        if oddPeak(0,peaks):
605            A[0] *=2
606        if oddPeak(1,peaks):
607            A[1] *=2
608        if oddPeak(2,peaks):
609            A[2] *=2
610    return A
611   
612def getDmin(peaks):
613    'needs a doc string'
614    return peaks[-1][-2]
615   
616def getDmax(peaks):
617    'needs a doc string'
618    return peaks[0][-2]
619   
620def refinePeaksZ(peaks,wave,ibrav,A,Zero,ZeroRef):
621    'needs a doc string'
622    dmin = getDmin(peaks)
623    OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,ZeroRef)
624    Peaks = np.array(peaks).T
625    H = Peaks[4:7]
626    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave))
627    peaks = Peaks.T   
628    HKL = G2lat.GenHBravais(dmin,ibrav,A)
629    M20,X20 = calc_M20(peaks,HKL)
630    return len(HKL),M20,X20,Aref,Z
631   
632def refinePeaksZSS(peaks,wave,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef):
633    'needs a doc string'
634    dmin = getDmin(peaks)
635    OK,smin,Aref,Vref,Z,result = FitHKLZSS(wave,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef)
636    Peaks = np.array(peaks).T
637    H = Peaks[4:8]
638    Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqZSS(H,Aref,Vref,Z,Peaks[0],wave))  #H,A,vec,Z,tth,lam
639    peaks = Peaks.T   
640    HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref)
641    M20,X20 = calc_M20SS(peaks,HKL)
642    return len(HKL),M20,X20,Aref,Vref,Z
643   
644def refinePeaksT(peaks,difC,ibrav,A,Zero,ZeroRef):
645    'needs a doc string'
646    dmin = getDmin(peaks)
647    OK,smin,Aref,Z,result = FitHKLT(difC,ibrav,peaks,A,Zero,ZeroRef)
648    Peaks = np.array(peaks).T
649    H = Peaks[4:7]
650    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqT(H,Aref,Z,Peaks[0],difC))
651    peaks = Peaks.T   
652    HKL = G2lat.GenHBravais(dmin,ibrav,A)
653    M20,X20 = calc_M20(peaks,HKL)
654    return len(HKL),M20,X20,Aref,Z
655   
656def refinePeaks(peaks,ibrav,A):
657    'needs a doc string'
658    dmin = getDmin(peaks)
659    smin = 1.0e10
660    pwr = 8
661    maxTries = 10
662    OK = False
663    tries = 0
664    HKL = G2lat.GenHBravais(dmin,ibrav,A)
665    while len(HKL) > 2 and IndexPeaks(peaks,HKL)[0]:
666        Pwr = pwr - (tries % 2)
667        HKL = []
668        tries += 1
669        osmin = smin
670        oldA = A[:]
671        Vold = G2lat.calc_V(oldA)
672        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
673        Vnew = G2lat.calc_V(A)
674        if Vnew > 2.0*Vold or Vnew < 2.:
675            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
676            OK = False
677            continue
678        try:
679            HKL = G2lat.GenHBravais(dmin,ibrav,A)
680        except FloatingPointError:
681            A = oldA
682            OK = False
683            break
684        if len(HKL) == 0: break                         #absurd cell obtained!
685        rat = (osmin-smin)/smin
686        if abs(rat) < 1.0e-5 or not OK: break
687        if tries > maxTries: break
688    if OK:
689        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
690        Peaks = np.array(peaks).T
691        H = Peaks[4:7]
692        try:
693            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
694            peaks = Peaks.T
695        except FloatingPointError:
696            A = oldA
697       
698    M20,X20 = calc_M20(peaks,HKL)
699    return len(HKL),M20,X20,A
700       
701def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
702    'needs a doc string'
703# dlg & ncMax are used for wx progress bar
704# A != 0 find the best A near input A,
705# A = 0 for random cell, volume normalized to V1;
706# returns number of generated hkls, M20, X20 & A for best found
707    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
708    dmin = getDmin(peaks)-0.05
709    amin = 2.5
710    amax = 5.*getDmax(peaks)
711    Asave = []
712    GoOn = True
713    if A:
714        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
715        if len(HKL) > mHKL[ibrav]:
716            peaks = IndexPeaks(peaks,HKL)[1]
717            Asave.append([calc_M20(peaks,HKL),A[:]])
718    tries = 0
719    while tries < Ntries:
720        if A:
721            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
722            if ibrav in [11,12,13]:         #monoclinic & triclinic
723                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
724        else:
725            Abeg = ranAbyV(ibrav,amin,amax,V1)
726        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
727       
728        if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]:
729            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
730            Asave.append([calc_M20(peaks,HKL),Aref[:]])
731            if ibrav == 9:                          #C-centered orthorhombic
732                for i in range(2):
733                    Abeg = rotOrthoA(Abeg[:])
734                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
735                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
736                    peaks = IndexPeaks(peaks,HKL)[1]
737                    Asave.append([calc_M20(peaks,HKL),Aref[:]])
738            elif ibrav == 11:                      #C-centered monoclinic
739                Abeg = swapMonoA(Abeg[:])
740                Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
741                HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
742                peaks = IndexPeaks(peaks,HKL)[1]
743                Asave.append([calc_M20(peaks,HKL),Aref[:]])
744        else:
745            break
746        Nc = len(HKL)
747        if Nc >= ncMax:
748            GoOn = False
749        elif dlg:
750            GoOn = dlg.Update(Nc)[0]
751            if not GoOn:
752                break
753        tries += 1   
754    X = sortM20(Asave)
755    if X:
756        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
757        return GoOn,Lhkl,M20,X20,A
758       
759    else:
760        return GoOn,0,0,0,0
761       
762def monoCellReduce(ibrav,A):
763    'needs a doc string'
764    a,b,c,alp,bet,gam = G2lat.A2cell(A)
765    G,g = G2lat.A2Gmat(A)
766    if ibrav in [11]:
767        u = [0,0,-1]
768        v = [1,0,2]
769        anew = math.sqrt(np.dot(np.dot(v,g),v))
770        if anew < a:
771            cang = np.dot(np.dot(u,g),v)/(anew*c)
772            beta = acosd(-abs(cang))
773            A = G2lat.cell2A([anew,b,c,90,beta,90])
774    else:
775        u = [-1,0,0]
776        v = [1,0,1]
777        cnew = math.sqrt(np.dot(np.dot(v,g),v))
778        if cnew < c:
779            cang = np.dot(np.dot(u,g),v)/(a*cnew)
780            beta = acosd(-abs(cang))
781            A = G2lat.cell2A([a,b,cnew,90,beta,90])
782    return A
783
784def DoIndexPeaks(peaks,controls,bravais):
785    'needs a doc string'
786   
787    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
788    amin = 2.5
789    amax = 5.0*getDmax(peaks)
790    dmin = getDmin(peaks)-delt
791    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
792        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
793        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
794    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
795    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
796    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
797    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
798    Nobs = len(peaks)
799    zero,ncno = controls[1:3]
800    ncMax = Nobs*ncno
801    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
802    print "%s %.4f %s %d %s %d" % ('Zero =',zero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
803    cells = []
804    for ibrav in range(14):
805        begin = time.time()
806        if bravais[ibrav]:
807            print 'cell search for ',bravaisNames[ibrav]
808            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
809            V1 = controls[3]
810            bestM20 = 0
811            topM20 = 0
812            cycle = 0
813            while cycle < 5:
814                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, 
815                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
816                screenSize = wx.ClientDisplayRect()
817                Size = dlg.GetSize()
818                dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
819                try:
820                    GoOn = True
821                    while GoOn:                                                 #Loop over increment of volume
822                        N2 = 0
823                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
824                            if ibrav > 2:
825                                if not N2:
826                                    A = []
827                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
828                                if A:
829                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
830                            else:
831                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
832                            if Nc >= ncMax:
833                                GoOn = False
834                                break
835                            elif 3*Nc < Nobs:
836                                N2 = 10
837                                break
838                            else:
839                                if not GoOn:
840                                    break
841                                if M20 > 1.0:
842                                    bestM20 = max(bestM20,M20)
843                                    A = halfCell(ibrav,A[:],peaks)
844                                    if ibrav in [12]:
845                                        A = monoCellReduce(ibrav,A[:])
846                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
847                                    peaks = IndexPeaks(peaks,HKL)[1]
848                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
849                                    V = G2lat.calc_V(A)
850                                    print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
851                                    if M20 >= 2.0:
852                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False])
853                            if not GoOn:
854                                break
855                            N2 += 1
856                        if ibrav < 11:
857                            V1 *= 1.1
858                        elif ibrav in range(11,14):
859                            V1 *= 1.05
860                        if not GoOn:
861                            if bestM20 > topM20:
862                                topM20 = bestM20
863                                if cells:
864                                    V1 = cells[0][9]
865                                else:
866                                    V1 = 25
867                                ncMax += Nobs
868                                cycle += 1
869                                print 'Restart search, new Max Nc = ',ncMax
870                            else:
871                                cycle = 10
872                finally:
873                    dlg.Destroy()
874            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
875                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin))
876           
877    if cells:
878        return True,dmin,cells
879    else:
880        return False,0,[]
881       
882       
883NeedTestData = True
884def TestData():
885    'needs a doc string'
886    array = np.array
887    global NeedTestData
888    NeedTestData = False
889    global TestData
890    TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
891        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], 
892        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], 
893        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], 
894        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], 
895        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], 
896        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], 
897        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
898        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], 
899        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], 
900        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], 
901        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], 
902        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], 
903        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], 
904        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
905        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], 
906        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], 
907        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], 
908        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], 
909        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], 
910        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], 
911        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], 
912        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], 
913        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], 
914        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], 
915        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], 
916        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], 
917        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], 
918        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
919        ]
920    global TestData2
921    TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
922        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], 
923        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], 
924        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], 
925        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], 
926        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], 
927        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], 
928        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], 
929        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], 
930        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], 
931        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], 
932        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
933        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], 
934        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], 
935        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], 
936        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], 
937        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], 
938        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], 
939        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], 
940        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], 
941        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], 
942        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], 
943        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], 
944        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], 
945        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], 
946        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], 
947        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], 
948        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], 
949        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
950        ]
951#
952def test0():
953    if NeedTestData: TestData()
954    msg = 'test FitHKL'
955    ibrav,cell,bestcell,Pwr,peaks = TestData
956    print 'best cell:',bestcell
957    print 'old cell:',cell
958    Peaks = np.array(peaks)
959    HKL = Peaks[4:7]
960    print calc_M20(peaks,HKL)
961    A = G2lat.cell2A(cell)
962    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
963    print 'new cell:',G2lat.A2cell(A)   
964    print 'x:',result[0]
965    print 'cov_x:',result[1]
966    print 'infodict:'
967    for item in result[2]:
968        print item,result[2][item]
969    print 'msg:',result[3]
970    print 'ier:',result[4]
971    result = refinePeaks(peaks,ibrav,A)
972    N,M20,X20,A = result
973    print 'refinePeaks:',N,M20,X20,G2lat.A2cell(A)
974    print 'compare bestcell:',bestcell
975#
976def test1():
977    if NeedTestData: TestData()
978    msg = 'test FitHKL'
979    ibrav,A,Pwr,peaks = TestData2
980    print 'bad cell:',G2lat.A2cell(A)
981    print 'FitHKL'
982    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
983    result = refinePeaks(peaks,ibrav,A)
984    N,M20,X20,A = result
985    print 'refinePeaks:',N,M20,X20,A
986#    Peaks = np.array(peaks)
987#    HKL = Peaks[4:7]
988#    print calc_M20(peaks,HKL)
989#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
990#    print 'new cell:',G2lat.A2cell(A)   
991#    print 'x:',result[0]
992#    print 'cov_x:',result[1]
993#    print 'infodict:'
994#    for item in result[2]:
995#        print item,result[2][item]
996#    print 'msg:',result[3]
997#    print 'ier:',result[4]
998   
999#
1000if __name__ == '__main__':
1001    test0()
1002    test1()
1003#    test2()
1004#    test3()
1005#    test4()
1006#    test5()
1007#    test6()
1008#    test7()
1009#    test8()
1010    print "OK"
Note: See TracBrowser for help on using the repository browser.