source: trunk/GSASIIpeak.py @ 243

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

GSASIIpeak.py - added factorize for future use in finding good numbers for fft
GSASIIimgGUI.py - minor tweak to q-value formatting & make small angle name 'SASD'
GSASIIIO.py - major refactor of GetTifData?

  • Property svn:keywords set to Date Author Revision URL Id
File size: 11.2 KB
Line 
1#GSASII peak fitting module
2########### SVN repository information ###################
3# $Date: 2011-02-08 21:47:07 +0000 (Tue, 08 Feb 2011) $
4# $Author: vondreele $
5# $Revision: 243 $
6# $URL: trunk/GSASIIpeak.py $
7# $Id: GSASIIpeak.py 243 2011-02-08 21:47:07Z vondreele $
8########### SVN repository information ###################
9import sys
10import math
11import wx
12import time
13import numpy as np
14import numpy.linalg as nl
15import GSASIIpath
16import pypowder as pyp              #assumes path has been amended to include correctr bin directory
17import GSASIIplot as G2plt
18
19# trig functions in degrees
20sind = lambda x: math.sin(x*math.pi/180.)
21asind = lambda x: 180.*math.asin(x)/math.pi
22tand = lambda x: math.tan(x*math.pi/180.)
23atand = lambda x: 180.*math.atan(x)/math.pi
24atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
25cosd = lambda x: math.cos(x*math.pi/180.)
26acosd = lambda x: 180.*math.acos(x)/math.pi
27rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
28#numpy versions
29npsind = lambda x: np.sin(x*np.pi/180.)
30npasind = lambda x: 180.*np.arcsin(x)/math.pi
31npcosd = lambda x: np.cos(x*math.pi/180.)
32nptand = lambda x: np.tan(x*math.pi/180.)
33npatand = lambda x: 180.*np.arctan(x)/np.pi
34npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35
36def factorize(num):
37    ''' Provide prime number factors for integer num
38    Returns dictionary of prime factors (keys) & power for each (data)
39    '''
40    factors = {}
41    orig = num
42
43    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
44    i, sqi = 2, 4
45    while sqi <= num:
46        while not num%i:
47            num /= i
48            factors[i] = factors.get(i, 0) + 1
49
50        sqi += 2*i + 1
51        i += 1
52
53    if num != 1 and num != orig:
54        factors[num] = factors.get(num, 0) + 1
55
56    if factors:
57        return factors
58    else:
59        return {num:1}          #a prime number!
60           
61def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
62    ''' Provide list of optimal data sizes for FFT calculations
63    Input:
64        nmin: minimum data size >= 1
65        nmax: maximum data size > nmin
66        thresh: maximum prime factor allowed
67    Returns:
68        list of data sizes where the maximum prime factor is < thresh
69    ''' 
70    plist = []
71    nmin = max(1,nmin)
72    nmax = max(nmin+1,nmax)
73    for p in range(nmin,nmax):
74        if max(factorize(p).keys()) < thresh:
75            plist.append(p)
76    return plist
77
78def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
79    # returns value(esd) string; nTZ=True for no trailing zeros
80    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
81    #get the 2 significant digits in the esd
82    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
83    #get the number of digits to represent them
84    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
85   
86    mdec = lambda esd: -int(math.log10(abs(esd)))
87    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
88    if esd > 0:
89        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
90        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
91        return fmt%(value,int(esd*10**(mdec(esd)+1)))
92    elif esd < 0:
93         return str(round(value,mdec(esd)))
94    else:
95        text = "%F"%(value)
96        if nTZ:
97            return text.rstrip('0')
98        else:
99            return text
100
101       
102#GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model       
103
104def DoPeakFit(peaks,background,limits,inst,data):
105   
106    def backgroundPrint(background,sigback):
107        if background[1]:
108            print 'Background coefficients for',background[0],'function'
109            ptfmt = "%12.5f"
110            ptstr =  'values:'
111            sigstr = 'esds  :'
112            for i,back in enumerate(background[3:]):
113                ptstr += ptfmt % (back)
114                sigstr += ptfmt % (sigback[i+3])
115            print ptstr
116            print sigstr
117        else:
118            print 'Background not refined'
119   
120    def instPrint(instVal,siginst,insLabels):
121        print 'Instrument Parameters:'
122        ptfmt = "%12.6f"
123        ptlbls = 'names :'
124        ptstr =  'values:'
125        sigstr = 'esds  :'
126        for i,value in enumerate(instVal):
127            ptlbls += "%s" % (insLabels[i].center(12))
128            ptstr += ptfmt % (value)
129            if siginst[i]:
130                sigstr += ptfmt % (siginst[i])
131            else:
132                sigstr += 12*' '
133        print ptlbls
134        print ptstr
135        print sigstr
136   
137    def peaksPrint(peaks,sigpeaks):
138        print 'Peak list='
139
140    begin = time.time()
141    Np = len(peaks[0])
142    DataType = inst[1][0]
143    instVal = inst[1][1:]
144    Insref = inst[2][1:]
145    insLabels = inst[3][1:]
146    Ka2 = False
147    Ioff = 3
148    if len(instVal) == 12:
149        lamratio = instVal[1]/instVal[0]
150        Ka2 = True
151        Ioff = 5
152    insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
153    for peak in peaks:
154        dip = []
155        dip.append(tand(peak[0]/2.0)**2)
156        dip.append(tand(peak[0]/2.0))
157        dip.append(1.0/cosd(peak[0]/2.0))
158        dip.append(tand(peak[0]/2.0))
159        peak.append(dip)
160    B = background[2]
161    bcof = background[3:3+B]
162    Bv = 0
163    if background[1]:
164        Bv = B
165    x,y,w,yc,yb,yd = data               #these are numpy arrays!
166    V = []
167    A = []
168    swobs = 0.0
169    smin = 0.0
170    first = True
171    GoOn = True
172    Go = True
173    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
174        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
175    screenSize = wx.DisplaySize()
176    Size = dlg.GetSize()
177    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
178    try:
179        i = 0
180        for xi in x :
181            Go = dlg.Update(i)[0]
182            if GoOn:
183                GoOn = Go
184            if limits[0] <= xi <= limits[1]:
185                yb[i] = 0.0
186                dp = []
187                for j in range(B):
188                    t = (xi-limits[0])**j
189                    yb[i] += t*bcof[j]
190                    if background[1]:
191                        dp.append(t)
192                yc[i] = yb[i]
193                Iv = 0
194                for j in range(6):
195                    if insref[j]:
196                        dp.append(0.0)
197                        Iv += 1
198                for peak in peaks:
199                    dip = peak[-1]
200                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0)
201                    yc[i] += f[0]*peak[2]
202                    if f[0] > 0.0:
203                        j = 0
204                        if insref[0]:              #U
205                            dp[Bv+j] += f[3]*dip[0]
206                            j += 1
207                        if insref[1]:              #V
208                            dp[Bv+j] += f[3]*dip[1]
209                            j += 1
210                        if insref[2]:              #W
211                            dp[Bv+j] += f[3]
212                            j += 1
213                        if insref[3]:              #X
214                            dp[Bv+j] += f[4]*dip[2]
215                            j += 1
216                        if insref[4]:              #Y
217                            dp[Bv+j] += f[4]*dip[3]
218                            j += 1
219                        if insref[5]:              #SH/L
220                            dp[Bv+j] += f[5]
221                    if Ka2:
222                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
223                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0)
224                       yc[i] += f2[0]*peak[2]*instVal[3]
225                       if f[0] > 0.0:
226                           j = 0
227                           if insref[0]:              #U
228                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
229                               j += 1
230                           if insref[1]:              #V
231                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
232                               j += 1
233                           if insref[2]:              #W
234                               dp[Bv+j] += f2[3]*instVal[3]
235                               j += 1
236                           if insref[3]:              #X
237                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
238                               j += 1
239                           if insref[4]:              #Y
240                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
241                               j += 1
242                           if insref[5]:              #SH/L
243                               dp[Bv+j] += f2[5]*instVal[3]                       
244                    for j in range(0,Np,2):
245                        if peak[j+1]: dp.append(f[j/2+1])
246                yd[i] = y[i]-yc[i]
247                swobs += w[i]*y[i]**2
248                t2 = w[i]*yd[i]
249                smin += t2*yd[i]
250                if first:
251                    first = False
252                    M = len(dp)
253                    A = np.zeros(shape=(M,M))
254                    V = np.zeros(shape=(M))
255                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
256            i += 1
257    finally:
258        dlg.Destroy()
259    Rwp = smin/swobs
260    Rwp = math.sqrt(Rwp)*100.0
261    norm = np.diag(A)
262    for i,elm in enumerate(norm):
263        if elm <= 0.0:
264            print norm
265            return False,0,0,0,False
266    for i in xrange(len(V)):
267        norm[i] = 1.0/math.sqrt(norm[i])
268        V[i] *= norm[i]
269        a = A[i]
270        for j in xrange(len(V)):
271            a[j] *= norm[i]
272    A = np.transpose(A)
273    for i in xrange(len(V)):
274        a = A[i]
275        for j in xrange(len(V)):
276            a[j] *= norm[i]
277    b = nl.solve(A,V)
278    A = nl.inv(A)
279    sig = np.diag(A)
280    for i in xrange(len(V)):
281        b[i] *= norm[i]
282        sig[i] *= norm[i]*norm[i]
283        sig[i] = math.sqrt(abs(sig[i]))
284    sigback = [0,0,0]
285    for j in range(Bv):
286        background[j+3] += b[j]
287        sigback.append(sig[j])
288    backgroundPrint(background,sigback)
289    k = 0
290    delt = []
291    if Ka2:
292        siginst = [0,0,0,0,0]
293    else:
294        siginst = [0,0,0]
295    for j in range(6):
296        if insref[j]:
297            instVal[j+Ioff] += b[Bv+k]
298            siginst.append(sig[Bv+k])
299            delt.append(b[Bv+k])
300            k += 1
301        else:
302            delt.append(0.0)
303            siginst.append(0.0)
304    delt.append(0.0)                    #dummies for azm
305    siginst.append(0.0)
306    instPrint(instVal,siginst,insLabels)
307    inst[1] = [DataType,]
308    for val in instVal:
309        inst[1].append(val)
310    B = Bv+Iv
311    for peak in peaks:
312        del peak[-1]                        # remove dip from end
313        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
314        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
315        for j in range(0,len(peak[:-1]),2):
316            if peak[j+1]: 
317                peak[j] += b[B]
318                B += 1
319        peak[4] += delsig
320        if peak[4] < 0.0:
321            print 'ERROR - negative sigma'
322            return False,0,0,0,False           
323        peak[6] += delgam
324        if peak[6] < 0.0:
325            print 'ERROR - negative gamma'
326            return False,0,0,0,False
327    runtime = time.time()-begin   
328    data = [x,y,w,yc,yb,yd]
329    return True,smin,Rwp,runtime,GoOn
330   
331           
332       
Note: See TracBrowser for help on using the repository browser.