source: trunk/GSASIIindex.py @ 1572

Last change on this file since 1572 was 1572, checked in by vondreele, 9 years ago

complete SS indexing, apply hklm extinction rules
cleanup indexing, cell refine, load cell, make new phase, calibration, etc.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 29.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII cell indexing program: variation on that of A. Coehlo
3#   includes cell refinement from peak positions (not zero as yet)
4########### SVN repository information ###################
5# $Date: 2014-11-17 22:37:02 +0000 (Mon, 17 Nov 2014) $
6# $Author: vondreele $
7# $Revision: 1572 $
8# $URL: trunk/GSASIIindex.py $
9# $Id: GSASIIindex.py 1572 2014-11-17 22:37:02Z 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: 1572 $")
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,peaks
233    hklds = list(np.array(HKL).T[-2])+[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[-1] >= 0:                                 # peak already assigned - test if this one better
248                opeak = peaks[hkl[-1]]
249                dold = abs(opeak[7]-hkl[-2])
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[-1] = ipk
257            peak[4:7] = hkl[:3]
258            peak[8] = hkl[-2]                                # 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,peaks
268    else:
269        return False,peaks  #nothing indexed!
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 FitHKLT(difC,ibrav,peaks,A,Z,Zref,Pwr):
342    'needs a doc string'
343   
344    def errFit(values,ibrav,d,H,tof,difC,Zref,Pwr):
345        Zero = Z
346        if Zref:   
347            Zero = values[-1]
348        A = Values2A(ibrav,values[:6])
349        Qo = 1./d**2
350        Qc = G2lat.calc_rDsqT(H,A,Zero,tof,difC)
351        return (Qo-Qc)*d**Pwr
352   
353    Peaks = np.array(peaks).T
354   
355    values = A2values(ibrav,A)
356    if Zref:
357        values.append(Z)
358    result = so.leastsq(errFit,values,full_output=True,ftol=0.0001,factor=0.001,
359        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Zref,Pwr))
360    A = Values2A(ibrav,result[0][:6])
361    if Zref:
362        Z = result[0][-1]
363   
364    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Zref,Pwr)**2),A,Z,result
365               
366def rotOrthoA(A):
367    'needs a doc string'
368    return [A[1],A[2],A[0],0,0,0]
369   
370def swapMonoA(A):
371    'needs a doc string'
372    return [A[2],A[1],A[0],0,A[4],0]
373   
374def oddPeak(indx,peaks):
375    'needs a doc string'
376    noOdd = True
377    for peak in peaks:
378        H = peak[4:7]
379        if H[indx] % 2:
380            noOdd = False
381    return noOdd
382   
383def halfCell(ibrav,A,peaks):
384    'needs a doc string'
385    if ibrav in [0,1,2]:
386        if oddPeak(0,peaks):
387            A[0] *= 2
388            A[1] = A[2] = A[0]
389    elif ibrav in [3,4,5,6]:
390        if oddPeak(0,peaks):
391            A[0] *= 2
392            A[1] = A[0]
393        if oddPeak(2,peaks):
394            A[2] *=2
395    else:
396        if oddPeak(0,peaks):
397            A[0] *=2
398        if oddPeak(1,peaks):
399            A[1] *=2
400        if oddPeak(2,peaks):
401            A[2] *=2
402    return A
403   
404def getDmin(peaks):
405    'needs a doc string'
406    return peaks[-1][7]
407   
408def getDmax(peaks):
409    'needs a doc string'
410    return peaks[0][7]
411   
412def refinePeaksZ(peaks,wave,ibrav,A,Zero,ZeroRef):
413    'needs a doc string'
414    dmin = getDmin(peaks)
415    OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,ZeroRef,0)
416    Peaks = np.array(peaks).T
417    H = Peaks[4:7]
418    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave))
419    peaks = Peaks.T   
420    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
421    M20,X20 = calc_M20(peaks,HKL)
422    return len(HKL),M20,X20,Aref,Z
423   
424def refinePeaksT(peaks,difC,ibrav,A,Zero,ZeroRef):
425    'needs a doc string'
426    dmin = getDmin(peaks)
427    OK,smin,Aref,Z,result = FitHKLT(difC,ibrav,peaks,A,Zero,ZeroRef,0)
428    Peaks = np.array(peaks).T
429    H = Peaks[4:7]
430    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqT(H,Aref,Z,Peaks[0],difC))
431    peaks = Peaks.T   
432    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
433    M20,X20 = calc_M20(peaks,HKL)
434    return len(HKL),M20,X20,Aref,Z
435   
436def refinePeaks(peaks,ibrav,A):
437    'needs a doc string'
438    dmin = getDmin(peaks)
439    smin = 1.0e10
440    pwr = 8
441    maxTries = 10
442    OK = False
443    tries = 0
444    HKL = G2lat.GenHBravais(dmin,ibrav,A)
445    while len(HKL) > 2 and IndexPeaks(peaks,HKL)[0]:
446        Pwr = pwr - (tries % 2)
447        HKL = []
448        tries += 1
449        osmin = smin
450        oldA = A[:]
451        Vold = G2lat.calc_V(oldA)
452        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
453        Vnew = G2lat.calc_V(A)
454        if Vnew > 2.0*Vold or Vnew < 2.:
455            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
456            OK = False
457            continue
458        try:
459            HKL = G2lat.GenHBravais(dmin,ibrav,A)
460        except FloatingPointError:
461            A = oldA
462            OK = False
463            break
464        if len(HKL) == 0: break                         #absurd cell obtained!
465        rat = (osmin-smin)/smin
466        if abs(rat) < 1.0e-5 or not OK: break
467        if tries > maxTries: break
468    if OK:
469        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
470        Peaks = np.array(peaks).T
471        H = Peaks[4:7]
472        try:
473            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
474            peaks = Peaks.T
475        except FloatingPointError:
476            A = oldA
477       
478    M20,X20 = calc_M20(peaks,HKL)
479    return len(HKL),M20,X20,A
480       
481def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
482    'needs a doc string'
483# dlg & ncMax are used for wx progress bar
484# A != 0 find the best A near input A,
485# A = 0 for random cell, volume normalized to V1;
486# returns number of generated hkls, M20, X20 & A for best found
487    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
488    dmin = getDmin(peaks)-0.05
489    amin = 2.5
490    amax = 5.*getDmax(peaks)
491    Asave = []
492    GoOn = True
493    if A:
494        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
495        if len(HKL) > mHKL[ibrav]:
496            peaks = IndexPeaks(peaks,HKL)[1]
497            Asave.append([calc_M20(peaks,HKL),A[:]])
498    tries = 0
499    while tries < Ntries:
500        if A:
501            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
502            if ibrav in [11,12,13]:         #monoclinic & triclinic
503                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
504        else:
505            Abeg = ranAbyV(ibrav,amin,amax,V1)
506        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
507       
508        if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]:
509            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
510            Asave.append([calc_M20(peaks,HKL),Aref[:]])
511            if ibrav == 9:                          #C-centered orthorhombic
512                for i in range(2):
513                    Abeg = rotOrthoA(Abeg[:])
514                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
515                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
516                    peaks = IndexPeaks(peaks,HKL)[1]
517                    Asave.append([calc_M20(peaks,HKL),Aref[:]])
518            elif ibrav == 11:                      #C-centered monoclinic
519                Abeg = swapMonoA(Abeg[:])
520                Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
521                HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
522                peaks = IndexPeaks(peaks,HKL)[1]
523                Asave.append([calc_M20(peaks,HKL),Aref[:]])
524        else:
525            break
526        Nc = len(HKL)
527        if Nc >= ncMax:
528            GoOn = False
529        elif dlg:
530            GoOn = dlg.Update(Nc)[0]
531            if not GoOn:
532                break
533        tries += 1   
534    X = sortM20(Asave)
535    if X:
536        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
537        return GoOn,Lhkl,M20,X20,A
538       
539    else:
540        return GoOn,0,0,0,0
541       
542def monoCellReduce(ibrav,A):
543    'needs a doc string'
544    a,b,c,alp,bet,gam = G2lat.A2cell(A)
545    G,g = G2lat.A2Gmat(A)
546    if ibrav in [11]:
547        u = [0,0,-1]
548        v = [1,0,2]
549        anew = math.sqrt(np.dot(np.dot(v,g),v))
550        if anew < a:
551            cang = np.dot(np.dot(u,g),v)/(anew*c)
552            beta = acosd(-abs(cang))
553            A = G2lat.cell2A([anew,b,c,90,beta,90])
554    else:
555        u = [-1,0,0]
556        v = [1,0,1]
557        cnew = math.sqrt(np.dot(np.dot(v,g),v))
558        if cnew < c:
559            cang = np.dot(np.dot(u,g),v)/(a*cnew)
560            beta = acosd(-abs(cang))
561            A = G2lat.cell2A([a,b,cnew,90,beta,90])
562    return A
563
564def DoIndexPeaks(peaks,controls,bravais):
565    'needs a doc string'
566   
567    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
568    amin = 2.5
569    amax = 5.0*getDmax(peaks)
570    dmin = getDmin(peaks)-delt
571    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
572        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
573        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
574    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
575    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
576    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
577    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
578    Nobs = len(peaks)
579    zero,ncno = controls[1:3]
580    ncMax = Nobs*ncno
581    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
582    print "%s %.4f %s %d %s %d" % ('Zero =',zero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
583    cells = []
584    for ibrav in range(14):
585        begin = time.time()
586        if bravais[ibrav]:
587            print 'cell search for ',bravaisNames[ibrav]
588            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
589            V1 = controls[3]
590            bestM20 = 0
591            topM20 = 0
592            cycle = 0
593            while cycle < 5:
594                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, 
595                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
596                screenSize = wx.ClientDisplayRect()
597                Size = dlg.GetSize()
598                dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
599                try:
600                    GoOn = True
601                    while GoOn:                                                 #Loop over increment of volume
602                        N2 = 0
603                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
604                            if ibrav > 2:
605                                if not N2:
606                                    A = []
607                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
608                                if A:
609                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
610                            else:
611                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
612                            if Nc >= ncMax:
613                                GoOn = False
614                                break
615                            elif 3*Nc < Nobs:
616                                N2 = 10
617                                break
618                            else:
619                                if not GoOn:
620                                    break
621                                if M20 > 1.0:
622                                    bestM20 = max(bestM20,M20)
623                                    A = halfCell(ibrav,A[:],peaks)
624                                    if ibrav in [12]:
625                                        A = monoCellReduce(ibrav,A[:])
626                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
627                                    peaks = IndexPeaks(peaks,HKL)[1]
628                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
629                                    V = G2lat.calc_V(A)
630                                    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)
631                                    if M20 >= 2.0:
632                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False])
633                            if not GoOn:
634                                break
635                            N2 += 1
636                        if ibrav < 11:
637                            V1 *= 1.1
638                        elif ibrav in range(11,14):
639                            V1 *= 1.05
640                        if not GoOn:
641                            if bestM20 > topM20:
642                                topM20 = bestM20
643                                if cells:
644                                    V1 = cells[0][9]
645                                else:
646                                    V1 = 25
647                                ncMax += Nobs
648                                cycle += 1
649                                print 'Restart search, new Max Nc = ',ncMax
650                            else:
651                                cycle = 10
652                finally:
653                    dlg.Destroy()
654            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
655                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin))
656           
657    if cells:
658        return True,dmin,cells
659    else:
660        return False,0,[]
661       
662       
663NeedTestData = True
664def TestData():
665    'needs a doc string'
666    array = np.array
667    global NeedTestData
668    NeedTestData = False
669    global TestData
670    TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
671        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], 
672        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], 
673        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], 
674        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], 
675        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], 
676        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], 
677        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
678        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], 
679        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], 
680        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], 
681        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], 
682        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], 
683        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], 
684        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
685        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], 
686        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], 
687        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], 
688        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], 
689        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], 
690        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], 
691        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], 
692        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], 
693        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], 
694        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], 
695        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], 
696        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], 
697        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], 
698        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
699        ]
700    global TestData2
701    TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
702        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], 
703        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], 
704        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], 
705        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], 
706        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], 
707        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], 
708        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], 
709        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], 
710        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], 
711        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], 
712        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
713        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], 
714        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], 
715        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], 
716        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], 
717        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], 
718        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], 
719        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], 
720        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], 
721        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], 
722        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], 
723        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], 
724        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], 
725        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], 
726        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], 
727        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], 
728        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], 
729        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
730        ]
731#
732def test0():
733    if NeedTestData: TestData()
734    msg = 'test FitHKL'
735    ibrav,cell,bestcell,Pwr,peaks = TestData
736    print 'best cell:',bestcell
737    print 'old cell:',cell
738    Peaks = np.array(peaks)
739    HKL = Peaks[4:7]
740    print calc_M20(peaks,HKL)
741    A = G2lat.cell2A(cell)
742    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
743    print 'new cell:',G2lat.A2cell(A)   
744    print 'x:',result[0]
745    print 'cov_x:',result[1]
746    print 'infodict:'
747    for item in result[2]:
748        print item,result[2][item]
749    print 'msg:',result[3]
750    print 'ier:',result[4]
751    result = refinePeaks(peaks,ibrav,A)
752    N,M20,X20,A = result
753    print 'refinePeaks:',N,M20,X20,G2lat.A2cell(A)
754    print 'compare bestcell:',bestcell
755#
756def test1():
757    if NeedTestData: TestData()
758    msg = 'test FitHKL'
759    ibrav,A,Pwr,peaks = TestData2
760    print 'bad cell:',G2lat.A2cell(A)
761    print 'FitHKL'
762    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
763    result = refinePeaks(peaks,ibrav,A)
764    N,M20,X20,A = result
765    print 'refinePeaks:',N,M20,X20,A
766#    Peaks = np.array(peaks)
767#    HKL = Peaks[4:7]
768#    print calc_M20(peaks,HKL)
769#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
770#    print 'new cell:',G2lat.A2cell(A)   
771#    print 'x:',result[0]
772#    print 'cov_x:',result[1]
773#    print 'infodict:'
774#    for item in result[2]:
775#        print item,result[2][item]
776#    print 'msg:',result[3]
777#    print 'ier:',result[4]
778   
779#
780if __name__ == '__main__':
781    test0()
782    test1()
783#    test2()
784#    test3()
785#    test4()
786#    test5()
787#    test6()
788#    test7()
789#    test8()
790    print "OK"
Note: See TracBrowser for help on using the repository browser.