source: trunk/GSASIIindex.py @ 1585

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

remove Pwr from FitPeaksZ, SSZ & T - wasn't used
put in derivatives for FitPeaksZ, etc. - works better
reorder arrows & make the expand/shrink symmetric about view point
remove zero & refine? from SS Unit Cells display - not refinable with modulation coeff.

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