source: branch/MPbranch/GSASIIindex.py

Last change on this file was 429, checked in by vondreele, 10 years ago

set 'browser' as default help display
fix zero & zero offset to be consistent in peak fitting/indexing/cell refinement

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