source: trunk/GSASIIindex.py @ 1582

Last change on this file since 1582 was 1582, checked in by vondreele, 8 years ago

revamp cell refinement in Unit Celle List
fix peak indexing in plots & bug in sqrt plot
put in error trap if indexing is tried with no peaks

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