source: trunk/GSASIIindex.py @ 774

Last change on this file since 774 was 762, checked in by toby, 13 years ago

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