source: trunk/GSASIIpeak.py @ 264

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

begin implementation of pdf generation - work in progress
modify image to have azimuth=0 as vertical "up"
add textctrls for max and min image intensity
fix to some lattice parameter defaults for some Laue groups

  • Property svn:keywords set to Date Author Revision URL Id
File size: 15.2 KB
Line 
1#GSASII powder calculation module
2########### SVN repository information ###################
3# $Date: 2011-03-31 21:41:13 +0000 (Thu, 31 Mar 2011) $
4# $Author: vondreele $
5# $Revision: 264 $
6# $URL: trunk/GSASIIpeak.py $
7# $Id: GSASIIpeak.py 264 2011-03-31 21:41:13Z 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 Transmission(Geometry,Abs,Diam):
79#Calculate sample transmission
80#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
81#   Abs: absorption coeff in cm-1
82#   Diam: sample thickness/diameter in mm
83    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
84        MuR = Abs*Diam/5.0
85        if MuR <= 3.0:
86            T0 = 16/(3.*math.pi)
87            T1 = -0.045780
88            T2 = -0.02489
89            T3 = 0.003045
90            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
91            if T < -20.:
92                return 2.06e-9
93            else:
94                return math.exp(T)
95        else:
96            T1 = 1.433902
97            T2 = 0.013869+0.337894
98            T3 = 1.933433+1.163198
99            T4 = 0.044365-0.04259
100            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
101            return T/100.
102
103def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
104#Calculate sample absorption
105#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
106#   Abs: absorption coeff in cm-1
107#   Diam: sample thickness/diameter in mm
108#   Tth: 2-theta scattering angle - can be numpy array
109#   Phi: flat plate tilt angle - future
110#   Psi: flat plate tilt axis - future
111    MuR = Abs*Diam/5.0
112    Sth2 = npsind(Tth/2.0)**2
113    Cth2 = 1.-Sth2
114    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
115        if MuR < 3.0:
116            T0 = 16.0/(3*np.pi)
117            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
118                0.109561*np.sqrt(Sth2)-26.04556
119            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
120                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
121            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
122            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
123            return np.exp(Trns)
124        else:
125            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
126                10.02088*Sth2**3-3.36778*Sth2**4
127            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
128                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
129            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
130                0.13576*np.sqrt(Sth2)+1.163198
131            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
132            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
133            return Trns/100.
134    elif 'Bragg' in Geometry:
135        return 1.0
136    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
137        # and only defined for 2theta < 90
138        T1 = np.exp(-MuR)
139        T2 = np.exp(-MuR/(1.-2.*Sth2))
140        Tb = -2.*Abs*Sth2
141        return (T1-T2)/Tb
142    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
143        cth = npcosd(Tth/2.0)
144        return (Diam/cth)*np.exp(-MuR/cth)
145       
146def Polarization(Pola,Azm,Tth):
147#   Calculate x-ray polarization correction
148#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
149#   Azm: azimuthal angle e.g. 0.0 in plane of polarization(?)
150#   Tth: 2-theta scattering angle - can be numpy array
151    pass
152   
153       
154
155def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
156    # returns value(esd) string; nTZ=True for no trailing zeros
157    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
158    #get the 2 significant digits in the esd
159    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
160    #get the number of digits to represent them
161    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
162   
163    mdec = lambda esd: -int(math.log10(abs(esd)))
164    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
165    if esd > 0:
166        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
167        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
168        return fmt%(value,int(esd*10**(mdec(esd)+1)))
169    elif esd < 0:
170         return str(round(value,mdec(esd)))
171    else:
172        text = "%F"%(value)
173        if nTZ:
174            return text.rstrip('0')
175        else:
176            return text
177
178       
179#GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model       
180
181def DoPeakFit(peaks,background,limits,inst,data):
182   
183    def backgroundPrint(background,sigback):
184        if background[1]:
185            print 'Background coefficients for',background[0],'function'
186            ptfmt = "%12.5f"
187            ptstr =  'values:'
188            sigstr = 'esds  :'
189            for i,back in enumerate(background[3:]):
190                ptstr += ptfmt % (back)
191                sigstr += ptfmt % (sigback[i+3])
192            print ptstr
193            print sigstr
194        else:
195            print 'Background not refined'
196   
197    def instPrint(instVal,siginst,insLabels):
198        print 'Instrument Parameters:'
199        ptfmt = "%12.6f"
200        ptlbls = 'names :'
201        ptstr =  'values:'
202        sigstr = 'esds  :'
203        for i,value in enumerate(instVal):
204            ptlbls += "%s" % (insLabels[i].center(12))
205            ptstr += ptfmt % (value)
206            if siginst[i]:
207                sigstr += ptfmt % (siginst[i])
208            else:
209                sigstr += 12*' '
210        print ptlbls
211        print ptstr
212        print sigstr
213   
214    def peaksPrint(peaks,sigpeaks):
215        print 'Peak list='
216
217    begin = time.time()
218    Np = len(peaks[0])
219    DataType = inst[1][0]
220    instVal = inst[1][1:]
221    Insref = inst[2][1:]
222    insLabels = inst[3][1:]
223    Ka2 = False
224    Ioff = 3
225    if len(instVal) == 12:
226        lamratio = instVal[1]/instVal[0]
227        Ka2 = True
228        Ioff = 5
229    insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
230    for peak in peaks:
231        dip = []
232        dip.append(tand(peak[0]/2.0)**2)
233        dip.append(tand(peak[0]/2.0))
234        dip.append(1.0/cosd(peak[0]/2.0))
235        dip.append(tand(peak[0]/2.0))
236        peak.append(dip)
237    B = background[2]
238    bcof = background[3:3+B]
239    Bv = 0
240    if background[1]:
241        Bv = B
242    x,y,w,yc,yb,yd = data               #these are numpy arrays!
243    V = []
244    A = []
245    swobs = 0.0
246    smin = 0.0
247    first = True
248    GoOn = True
249    Go = True
250    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
251        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
252    screenSize = wx.DisplaySize()
253    Size = dlg.GetSize()
254    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
255    try:
256        i = 0
257        for xi in x :
258            Go = dlg.Update(i)[0]
259            if GoOn:
260                GoOn = Go
261            if limits[0] <= xi <= limits[1]:
262                yb[i] = 0.0
263                dp = []
264                for j in range(B):
265                    t = (xi-limits[0])**j
266                    yb[i] += t*bcof[j]
267                    if background[1]:
268                        dp.append(t)
269                yc[i] = yb[i]
270                Iv = 0
271                for j in range(6):
272                    if insref[j]:
273                        dp.append(0.0)
274                        Iv += 1
275                for peak in peaks:
276                    dip = peak[-1]
277                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0)
278                    yc[i] += f[0]*peak[2]
279                    if f[0] > 0.0:
280                        j = 0
281                        if insref[0]:              #U
282                            dp[Bv+j] += f[3]*dip[0]
283                            j += 1
284                        if insref[1]:              #V
285                            dp[Bv+j] += f[3]*dip[1]
286                            j += 1
287                        if insref[2]:              #W
288                            dp[Bv+j] += f[3]
289                            j += 1
290                        if insref[3]:              #X
291                            dp[Bv+j] += f[4]*dip[2]
292                            j += 1
293                        if insref[4]:              #Y
294                            dp[Bv+j] += f[4]*dip[3]
295                            j += 1
296                        if insref[5]:              #SH/L
297                            dp[Bv+j] += f[5]
298                    if Ka2:
299                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
300                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0)
301                       yc[i] += f2[0]*peak[2]*instVal[3]
302                       if f[0] > 0.0:
303                           j = 0
304                           if insref[0]:              #U
305                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
306                               j += 1
307                           if insref[1]:              #V
308                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
309                               j += 1
310                           if insref[2]:              #W
311                               dp[Bv+j] += f2[3]*instVal[3]
312                               j += 1
313                           if insref[3]:              #X
314                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
315                               j += 1
316                           if insref[4]:              #Y
317                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
318                               j += 1
319                           if insref[5]:              #SH/L
320                               dp[Bv+j] += f2[5]*instVal[3]                       
321                    for j in range(0,Np,2):
322                        if peak[j+1]: dp.append(f[j/2+1])
323                yd[i] = y[i]-yc[i]
324                swobs += w[i]*y[i]**2
325                t2 = w[i]*yd[i]
326                smin += t2*yd[i]
327                if first:
328                    first = False
329                    M = len(dp)
330                    A = np.zeros(shape=(M,M))
331                    V = np.zeros(shape=(M))
332                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
333            i += 1
334    finally:
335        dlg.Destroy()
336    Rwp = smin/swobs
337    Rwp = math.sqrt(Rwp)*100.0
338    norm = np.diag(A)
339    for i,elm in enumerate(norm):
340        if elm <= 0.0:
341            print norm
342            return False,0,0,0,False
343    for i in xrange(len(V)):
344        norm[i] = 1.0/math.sqrt(norm[i])
345        V[i] *= norm[i]
346        a = A[i]
347        for j in xrange(len(V)):
348            a[j] *= norm[i]
349    A = np.transpose(A)
350    for i in xrange(len(V)):
351        a = A[i]
352        for j in xrange(len(V)):
353            a[j] *= norm[i]
354    b = nl.solve(A,V)
355    A = nl.inv(A)
356    sig = np.diag(A)
357    for i in xrange(len(V)):
358        b[i] *= norm[i]
359        sig[i] *= norm[i]*norm[i]
360        sig[i] = math.sqrt(abs(sig[i]))
361    sigback = [0,0,0]
362    for j in range(Bv):
363        background[j+3] += b[j]
364        sigback.append(sig[j])
365    backgroundPrint(background,sigback)
366    k = 0
367    delt = []
368    if Ka2:
369        siginst = [0,0,0,0,0]
370    else:
371        siginst = [0,0,0]
372    for j in range(6):
373        if insref[j]:
374            instVal[j+Ioff] += b[Bv+k]
375            siginst.append(sig[Bv+k])
376            delt.append(b[Bv+k])
377            k += 1
378        else:
379            delt.append(0.0)
380            siginst.append(0.0)
381    delt.append(0.0)                    #dummies for azm
382    siginst.append(0.0)
383    instPrint(instVal,siginst,insLabels)
384    inst[1] = [DataType,]
385    for val in instVal:
386        inst[1].append(val)
387    B = Bv+Iv
388    for peak in peaks:
389        del peak[-1]                        # remove dip from end
390        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
391        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
392        for j in range(0,len(peak[:-1]),2):
393            if peak[j+1]: 
394                peak[j] += b[B]
395                B += 1
396        peak[4] += delsig
397        if peak[4] < 0.0:
398            print 'ERROR - negative sigma'
399            return False,0,0,0,False           
400        peak[6] += delgam
401        if peak[6] < 0.0:
402            print 'ERROR - negative gamma'
403            return False,0,0,0,False
404    runtime = time.time()-begin   
405    data = [x,y,w,yc,yb,yd]
406    return True,smin,Rwp,runtime,GoOn
407
408def ComputePDF(data,xydata):
409    for key in data:
410        print key,data[key]
411    #subtract backgrounds - if any
412    xydata['Sample corrected'] = xydata['Sample']
413    if 'Sample Bkg.' in xydata:
414        xydata['Sample corrected'][1][1] -= (xydata['Sample Bkg.'][1][1]+
415            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
416    if 'Container' in xydata:   
417        xydata['Sample corrected'][1][1] -= (xydata['Container'][1][1]+
418            data['Container']['Add'])*data['Container']['Mult']
419    if 'Container Bkg.' in xydata:
420        xydata['Sample corrected'][1][1] += (xydata['Container Bkg.'][1][1]+
421            data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
422   
423           
424       
425    return xydata['Sample corrected'],[]
426       
Note: See TracBrowser for help on using the repository browser.