source: trunk/GSASIIindex.py @ 342

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

major modifications to get 1st "working" version of Refine
Add GSASIImapvars.py
fix cell refinement in indexing window
single cycle for peak refinement

  • Property svn:keywords set to Date Author Revision URL Id
File size: 27.4 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-08-05 19:35:43 +0000 (Fri, 05 Aug 2011) $
5# $Author: vondreele $
6# $Revision: 342 $
7# $URL: trunk/GSASIIindex.py $
8# $Id: GSASIIindex.py 342 2011-08-05 19:35:43Z 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 IndexPeaks(peaks,HKL):
195    import bisect
196    N = len(HKL)
197    if N == 0: return False
198    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
199    hklds.sort()                                        # ascending sort - upper bound at end
200    hklmax = [0,0,0]
201    for ipk,peak in enumerate(peaks):
202        if peak[2]:
203            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
204            dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
205            dp = hklds[i]-peak[7]
206            pos = N-i                                       # reverse the order
207            if dp > dm: pos += 1                            # closer to upper than lower
208            hkl = HKL[pos]                                 # put in hkl
209            if hkl[4] >= 0:                                 # peak already assigned - test if this one better
210                opeak = peaks[hkl[4]]
211                dold = abs(opeak[7]-hkl[3])
212                dnew = min(dm,dp)
213                if dold > dnew:                             # new better - zero out old
214                    opeak[4:7] = [0,0,0]
215                    opeak[8] = 0.
216                else:                                       # old better - do nothing
217                    continue               
218            hkl[4] = ipk
219            peak[4:7] = hkl[:3]
220            peak[8] = hkl[3]                                # fill in d-calc
221    for peak in peaks:
222        peak[3] = False
223        if peak[2]:
224            if peak[8] > 0.:
225                for j in range(3):
226                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
227                peak[3] = True
228    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
229        return True
230    else:
231        return False
232       
233def Values2A(ibrav,values):
234    if ibrav in [0,1,2]:
235        return [values[0],values[0],values[0],0,0,0]
236    elif ibrav in [3,4]:
237        return [values[0],values[0],values[1],values[0],0,0]
238    elif ibrav in [5,6]:
239        return [values[0],values[0],values[1],0,0,0]
240    elif ibrav in [7,8,9,10]:
241        return [values[0],values[1],values[2],0,0,0]
242    elif ibrav in [11,12]:
243        return [values[0],values[1],values[2],0,values[3],0]
244    else:
245        return list(values)
246       
247def A2values(ibrav,A):
248    if ibrav in [0,1,2]:
249        return [A[0],]
250    elif ibrav in [3,4,5,6]:
251        return [A[0],A[2]]
252    elif ibrav in [7,8,9,10]:
253        return [A[0],A[1],A[2]]
254    elif ibrav in [11,12]:
255        return [A[0],A[1],A[2],A[4]]
256    else:
257        return A
258
259def FitHKL(ibrav,peaks,A,Pwr):
260               
261    def errFit(values,ibrav,d,H,Pwr):
262        A = Values2A(ibrav,values)
263        Qo = 1./d**2
264        Qc = G2lat.calc_rDsq(H,A)
265        return (Qo-Qc)*d**Pwr
266   
267    Peaks = np.array(peaks).T
268   
269    values = A2values(ibrav,A)
270    result = so.leastsq(errFit,values,full_output=True,ftol=0.0001,
271        args=(ibrav,Peaks[7],Peaks[4:7],Pwr))
272    A = Values2A(ibrav,result[0])
273    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
274           
275def FitHKLZ(wave,ibrav,peaks,A,Z,Pwr):
276   
277    def errFit(values,ibrav,d,H,tth,wave,Pwr):
278        A = Values2A(ibrav,values[:-1])
279        Z = values[-1]
280        Qo = 1./d**2
281        Qc = G2lat.calc_rDsqZ(H,A,Z,tth,wave)
282        return (Qo-Qc)*d**Pwr
283   
284    Peaks = np.array(peaks).T
285   
286    values = A2values(ibrav,A)
287    values.append(Z)
288    result = so.leastsq(errFit,values,full_output=True,ftol=0.001,
289        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr))
290    A = Values2A(ibrav,result[0][:-1])
291    Z = result[0][-1]
292   
293    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Pwr)**2),A,Z,result
294               
295def rotOrthoA(A):
296    return [A[1],A[2],A[0],0,0,0]
297   
298def swapMonoA(A):
299    return [A[2],A[1],A[0],0,A[4],0]
300   
301def oddPeak(indx,peaks):
302    noOdd = True
303    for peak in peaks:
304        H = peak[4:7]
305        if H[indx] % 2:
306            noOdd = False
307    return noOdd
308   
309def halfCell(ibrav,A,peaks):
310    if ibrav in [0,1,2]:
311        if oddPeak(0,peaks):
312            A[0] *= 2
313            A[1] = A[2] = A[0]
314    elif ibrav in [3,4,5,6]:
315        if oddPeak(0,peaks):
316            A[0] *= 2
317            A[1] = A[0]
318        if oddPeak(2,peaks):
319            A[2] *=2
320    else:
321        if oddPeak(0,peaks):
322            A[0] *=2
323        if oddPeak(1,peaks):
324            A[1] *=2
325        if oddPeak(2,peaks):
326            A[2] *=2
327    return A
328   
329def getDmin(peaks):
330    return peaks[-1][7]
331   
332def getDmax(peaks):
333    return peaks[0][7]
334   
335def refinePeaksZ(peaks,wave,ibrav,A,Zero):
336    dmin = getDmin(peaks)
337    OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,2)
338    Peaks = np.array(peaks).T
339    H = Peaks[4:7]
340    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave))
341    peaks = Peaks.T
342   
343    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
344    M20,X20 = calc_M20(peaks,HKL)
345    return len(HKL),M20,X20,Aref,Z
346   
347def refinePeaks(peaks,ibrav,A):
348    dmin = getDmin(peaks)
349    smin = 1.0e10
350    pwr = 8
351    maxTries = 10
352    OK = False
353    tries = 0
354    HKL = G2lat.GenHBravais(dmin,ibrav,A)
355    while len(HKL) > 2 and IndexPeaks(peaks,HKL):
356        Pwr = pwr - (tries % 2)
357        HKL = []
358        tries += 1
359        osmin = smin
360        oldA = A[:]
361        Vold = G2lat.calc_V(oldA)
362        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
363        Vnew = G2lat.calc_V(A)
364        if Vnew > 2.0*Vold or Vnew < 2.:
365            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
366            OK = False
367            continue
368        try:
369            HKL = G2lat.GenHBravais(dmin,ibrav,A)
370        except FloatingPointError:
371            A = oldA
372            OK = False
373            break
374        if len(HKL) == 0: break                         #absurd cell obtained!
375        rat = (osmin-smin)/smin
376        if abs(rat) < 1.0e-5 or not OK: break
377        if tries > maxTries: break
378    if OK:
379        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
380        Peaks = np.array(peaks).T
381        H = Peaks[4:7]
382        try:
383            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
384        except FloatingPointError:
385            print G2lat.calc_rDsq(H,A)
386            Peaks[8] = 1.0
387        peaks = Peaks.T
388       
389    M20,X20 = calc_M20(peaks,HKL)
390    return len(HKL),M20,X20,A
391       
392def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
393# dlg & ncMax are used for wx progress bar
394# A != 0 find the best A near input A,
395# A = 0 for random cell, volume normalized to V1;
396# returns number of generated hkls, M20, X20 & A for best found
397    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
398    dmin = getDmin(peaks)-0.05
399    amin = 2.5
400    amax = 5.*getDmax(peaks)
401    Asave = []
402    GoOn = True
403    if A:
404        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
405        if len(HKL) > mHKL[ibrav]:
406            IndexPeaks(peaks,HKL)
407            Asave.append([calc_M20(peaks,HKL),A[:]])
408    tries = 0
409    while tries < Ntries:
410        if A:
411            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
412            if ibrav in [11,12,13]:         #monoclinic & triclinic
413                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
414        else:
415            Abeg = ranAbyV(ibrav,amin,amax,V1)
416        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
417       
418        if IndexPeaks(peaks,HKL) and len(HKL) > mHKL[ibrav]:
419            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
420            Asave.append([calc_M20(peaks,HKL),Aref[:]])
421            if ibrav == 9:                          #C-centered orthorhombic
422                for i in range(2):
423                    Abeg = rotOrthoA(Abeg[:])
424                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
425                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
426                    IndexPeaks(peaks,HKL)
427                    Asave.append([calc_M20(peaks,HKL),Aref[:]])
428            elif ibrav == 11:                      #C-centered monoclinic
429                Abeg = swapMonoA(Abeg[:])
430                Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
431                HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
432                IndexPeaks(peaks,HKL)
433                Asave.append([calc_M20(peaks,HKL),Aref[:]])
434        else:
435            break
436        Nc = len(HKL)
437        if Nc >= ncMax:
438            GoOn = False
439        elif dlg:
440            GoOn = dlg.Update(Nc)[0]
441            if not GoOn:
442                break
443        tries += 1   
444    X = sortM20(Asave)
445    if X:
446        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
447        return GoOn,Lhkl,M20,X20,A
448       
449    else:
450        return GoOn,0,0,0,0
451       
452def monoCellReduce(ibrav,A):
453    a,b,c,alp,bet,gam = G2lat.A2cell(A)
454    G,g = G2lat.A2Gmat(A)
455    if ibrav in [11]:
456        u = [0,0,-1]
457        v = [1,0,2]
458        anew = math.sqrt(np.dot(np.dot(v,g),v))
459        if anew < a:
460            cang = np.dot(np.dot(u,g),v)/(anew*c)
461            beta = acosd(-abs(cang))
462            A = G2lat.cell2A([anew,b,c,90,beta,90])
463    else:
464        u = [-1,0,0]
465        v = [1,0,1]
466        cnew = math.sqrt(np.dot(np.dot(v,g),v))
467        if cnew < c:
468            cang = np.dot(np.dot(u,g),v)/(a*cnew)
469            beta = acosd(-abs(cang))
470            A = G2lat.cell2A([a,b,cnew,90,beta,90])
471    return A
472
473def DoIndexPeaks(peaks,inst,controls,bravais):
474   
475    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
476    amin = 2.5
477    amax = 5.0*getDmax(peaks)
478    dmin = getDmin(peaks)-delt
479    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
480        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
481        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
482    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
483    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
484    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
485    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
486    Nobs = len(peaks)
487    wave = inst[1]
488    if len(inst) > 10:
489        zero = inst[3]
490    else:
491        zero = inst[2]
492    print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
493    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
494    ifzero,maxzero,ncno = controls[:3]
495    ncMax = Nobs*ncno
496    print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
497    cells = []
498    for ibrav in range(14):
499        begin = time.time()
500        if bravais[ibrav]:
501            print 'cell search for ',bravaisNames[ibrav]
502            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
503            V1 = controls[3]
504            bestM20 = 0
505            topM20 = 0
506            cycle = 0
507            while cycle < 5:
508                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, 
509                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
510                screenSize = wx.ClientDisplayRect()
511                Size = dlg.GetSize()
512                dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
513                try:
514                    GoOn = True
515                    while GoOn:                                                 #Loop over increment of volume
516                        N2 = 0
517                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
518                            if ibrav > 2:
519                                if not N2:
520                                    A = []
521                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
522                                if A:
523                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
524                            else:
525                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
526                            if Nc >= ncMax:
527                                GoOn = False
528                                break
529                            elif 3*Nc < Nobs:
530                                N2 = 10
531                                break
532                            else:
533                                if not GoOn:
534                                    break
535                                if M20 > 1.0:
536                                    bestM20 = max(bestM20,M20)
537                                    A = halfCell(ibrav,A[:],peaks)
538                                    if ibrav in [12]:
539                                        A = monoCellReduce(ibrav,A[:])
540                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
541                                    IndexPeaks(peaks,HKL)
542                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
543                                    V = G2lat.calc_V(A)
544                                    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)
545                                    if M20 >= 2.0:
546                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
547                            if not GoOn:
548                                break
549                            N2 += 1
550                        if ibrav < 11:
551                            V1 *= 1.1
552                        elif ibrav in range(11,14):
553                            V1 *= 1.05
554                        if not GoOn:
555                            if bestM20 > topM20:
556                                topM20 = bestM20
557                                if cells:
558                                    V1 = cells[0][9]
559                                else:
560                                    V1 = 25
561                                ncMax += Nobs
562                                cycle += 1
563                                print 'Restart search, new Max Nc = ',ncMax
564                            else:
565                                cycle = 10
566                finally:
567                    dlg.Destroy()
568            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
569                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin))
570           
571    if cells:
572        cells = sortM20(cells)
573        cells[0][-1] = True
574        return True,dmin,cells
575    else:
576        return False,0,0
577       
578       
579NeedTestData = True
580def TestData():
581    array = np.array
582    global NeedTestData
583    NeedTestData = False
584    global TestData
585    TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
586        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], 
587        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], 
588        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], 
589        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], 
590        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], 
591        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], 
592        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
593        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], 
594        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], 
595        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], 
596        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], 
597        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], 
598        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], 
599        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
600        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], 
601        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], 
602        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], 
603        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], 
604        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], 
605        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], 
606        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], 
607        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], 
608        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], 
609        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], 
610        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], 
611        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], 
612        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], 
613        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
614        ]
615    global TestData2
616    TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
617        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], 
618        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], 
619        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], 
620        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], 
621        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], 
622        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], 
623        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], 
624        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], 
625        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], 
626        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], 
627        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
628        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], 
629        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], 
630        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], 
631        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], 
632        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], 
633        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], 
634        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], 
635        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], 
636        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], 
637        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], 
638        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], 
639        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], 
640        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], 
641        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], 
642        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], 
643        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], 
644        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
645        ]
646#
647def test0():
648    if NeedTestData: TestData()
649    msg = 'test FitHKL'
650    ibrav,cell,bestcell,Pwr,peaks = TestData
651    print 'best cell:',bestcell
652    print 'old cell:',cell
653    Peaks = np.array(peaks)
654    HKL = Peaks[4:7]
655    print calc_M20(peaks,HKL)
656    A = G2lat.cell2A(cell)
657    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
658    print 'new cell:',G2lat.A2cell(A)   
659    print 'x:',result[0]
660    print 'cov_x:',result[1]
661    print 'infodict:'
662    for item in result[2]:
663        print item,result[2][item]
664    print 'msg:',result[3]
665    print 'ier:',result[4]
666    result = refinePeaks(peaks,ibrav,A)
667    N,M20,X20,A = result
668    print 'refinePeaks:',N,M20,X20,G2lat.A2cell(A)
669    print 'compare bestcell:',bestcell
670#
671def test1():
672    if NeedTestData: TestData()
673    msg = 'test FitHKL'
674    ibrav,A,Pwr,peaks = TestData2
675    print 'bad cell:',G2lat.A2cell(A)
676    print 'FitHKL'
677    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
678    result = refinePeaks(peaks,ibrav,A)
679    N,M20,X20,A = result
680    print 'refinePeaks:',N,M20,X20,A
681#    Peaks = np.array(peaks)
682#    HKL = Peaks[4:7]
683#    print calc_M20(peaks,HKL)
684#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
685#    print 'new cell:',G2lat.A2cell(A)   
686#    print 'x:',result[0]
687#    print 'cov_x:',result[1]
688#    print 'infodict:'
689#    for item in result[2]:
690#        print item,result[2][item]
691#    print 'msg:',result[3]
692#    print 'ier:',result[4]
693   
694#
695if __name__ == '__main__':
696    test0()
697    test1()
698#    test2()
699#    test3()
700#    test4()
701#    test5()
702#    test6()
703#    test7()
704#    test8()
705    print "OK"
Note: See TracBrowser for help on using the repository browser.