source: trunk/GSASIIpeak.py @ 82

Last change on this file since 82 was 82, checked in by vondreel, 12 years ago

split powder GUI stuff out of G2gd & put in new GSASIIpwdGUI.py - all seems to work

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