source: trunk/GSASIIpeak.py @ 196

Last change on this file since 196 was 196, checked in by vondreele, 13 years ago

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