source: trunk/GSASIIindex.py @ 1007

Last change on this file since 1007 was 939, checked in by toby, 12 years ago

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