source: trunk/GSASIIindex.py @ 1578

Last change on this file since 1578 was 1578, checked in by vondreele, 8 years ago

add wx.Yield in RefineBox?
refactor indexing for normal & superlattice PeaksPrint?
remove splitter window in Unit Cells Data display - didn't work in newer wxPythons

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