# source:trunk/GSASIIpeak.py@265

Last change on this file since 265 was 265, checked in by vondreele, 12 years ago

further progress on implementing pdf calculations
optionally put legends on the pdf plots
attempt implementation of a rotation of the azimuth ranges for multiazimuth integrations -not fully successful

• Property svn:keywords set to `Date Author Revision URL Id`
File size: 27.9 KB
Line
1#GSASII powder calculation module
2########### SVN repository information ###################
3# \$Date: 2011-04-19 20:12:37 +0000 (Tue, 19 Apr 2011) \$
4# \$Author: vondreele \$
5# \$Revision: 265 \$
6# \$URL: trunk/GSASIIpeak.py \$
7# \$Id: GSASIIpeak.py 265 2011-04-19 20:12:37Z vondreele \$
8########### SVN repository information ###################
9import sys
10import math
11import wx
12import time
13import numpy as np
14import scipy as sp
15import numpy.linalg as nl
16import scipy.interpolate as si
17import GSASIIpath
18import pypowder as pyp              #assumes path has been amended to include correctr bin directory
19import GSASIIplot as G2plt
20import GSASIIlattice as G2lat
21import GSASIIElem as G2elem
22import GSASIIgrid as G2gd
23
24# trig functions in degrees
25sind = lambda x: math.sin(x*math.pi/180.)
26asind = lambda x: 180.*math.asin(x)/math.pi
27tand = lambda x: math.tan(x*math.pi/180.)
28atand = lambda x: 180.*math.atan(x)/math.pi
29atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
30cosd = lambda x: math.cos(x*math.pi/180.)
31acosd = lambda x: 180.*math.acos(x)/math.pi
32rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
33#numpy versions
34npsind = lambda x: np.sin(x*np.pi/180.)
35npasind = lambda x: 180.*np.arcsin(x)/math.pi
36npcosd = lambda x: np.cos(x*math.pi/180.)
37nptand = lambda x: np.tan(x*math.pi/180.)
38npatand = lambda x: 180.*np.arctan(x)/np.pi
39npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
40npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
41npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
42
43def factorize(num):
44    ''' Provide prime number factors for integer num
45    Returns dictionary of prime factors (keys) & power for each (data)
46    '''
47    factors = {}
48    orig = num
49
50    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
51    i, sqi = 2, 4
52    while sqi <= num:
53        while not num%i:
54            num /= i
55            factors[i] = factors.get(i, 0) + 1
56
57        sqi += 2*i + 1
58        i += 1
59
60    if num != 1 and num != orig:
61        factors[num] = factors.get(num, 0) + 1
62
63    if factors:
64        return factors
65    else:
66        return {num:1}          #a prime number!
67
68def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
69    ''' Provide list of optimal data sizes for FFT calculations
70    Input:
71        nmin: minimum data size >= 1
72        nmax: maximum data size > nmin
73        thresh: maximum prime factor allowed
74    Returns:
75        list of data sizes where the maximum prime factor is < thresh
76    '''
77    plist = []
78    nmin = max(1,nmin)
79    nmax = max(nmin+1,nmax)
80    for p in range(nmin,nmax):
81        if max(factorize(p).keys()) < thresh:
82            plist.append(p)
83    return plist
84
85def Transmission(Geometry,Abs,Diam):
86#Calculate sample transmission
87#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
88#   Abs: absorption coeff in cm-1
89#   Diam: sample thickness/diameter in mm
90    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
91        MuR = Abs*Diam/20.0
92        if MuR <= 3.0:
93            T0 = 16/(3.*math.pi)
94            T1 = -0.045780
95            T2 = -0.02489
96            T3 = 0.003045
97            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
98            if T < -20.:
99                return 2.06e-9
100            else:
101                return math.exp(T)
102        else:
103            T1 = 1.433902
104            T2 = 0.013869+0.337894
105            T3 = 1.933433+1.163198
106            T4 = 0.044365-0.04259
107            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
108            return T/100.
109    elif 'plate' in Geometry:
110        MuR = Abs*Diam/10.
111        return math.exp(-MuR)
112    elif 'Bragg' in Geometry:
113        return 0.0
114
115def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
116#Calculate sample absorption
117#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
118#   Abs: absorption coeff in cm-1
119#   Diam: sample thickness/diameter in mm
120#   Tth: 2-theta scattering angle - can be numpy array
121#   Phi: flat plate tilt angle - future
122#   Psi: flat plate tilt axis - future
123    Sth2 = npsind(Tth/2.0)**2
124    Cth2 = 1.-Sth2
125    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
126        MuR = Abs*Diam/20.0
127        if MuR < 3.0:
128            T0 = 16.0/(3*np.pi)
129            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
130                0.109561*np.sqrt(Sth2)-26.04556
131            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
132                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
133            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
134            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
135            return np.exp(Trns)
136        else:
137            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
138                10.02088*Sth2**3-3.36778*Sth2**4
139            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
140                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
141            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
142                0.13576*np.sqrt(Sth2)+1.163198
143            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
144            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
145            return Trns/100.
146    elif 'Bragg' in Geometry:
147        return 1.0
148    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
149        # and only defined for 2theta < 90
150        MuR = Abs*Diam/10.0
151        T1 = np.exp(-MuR)
152        T2 = np.exp(-MuR/npcosd(Tth))
153        Tb = MuR-MuR/npcosd(Tth)
154        return (T2-T1)/Tb
155    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
156        MuR = Abs*Diam/10.0
157        cth = npcosd(Tth/2.0)
158        return np.exp(-MuR/cth)/cth
159
160def Polarization(Pola,Tth,Azm=0.0):
161#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
162#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
163#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
164#   Tth: 2-theta scattering angle - can be numpy array
165#       which (if either) of these is "right"?
166#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
167#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
168#    return (Pola*npcosd(Azm)**2+npsind(Azm)**2)*npcosd(Tth)**2+   \
169#        Pola*npsind(Azm)**2+npcosd(Azm)**2
170    return Pola*npcosd(Tth)**2+1.0
171
172def Oblique(ObCoeff,Tth):
173    if ObCoeff:
174        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
175    else:
176        return 1.0
177
178def Ruland(RulCoff,wave,Q,Compton):
179    C = 2.9978e8
180    D = 1.5e-3
181    hmc = 0.024262734687
182    sinth2 = (Q*wave/(4.0*np.pi))**2
183    dlam = (wave**2)*Compton*Q/C
184    dlam_c = 2.0*hmc*sinth2-D*wave**2
185    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
186
187def lambdaCompton(DetType,wave,Q):
188    hmc = 0.024262734687
189    if 'Image' in DetType:
190        return wave*(1.0+2.0*hmc*wave*(Q/(4.0*np.pi))**2)**3
191    else:
192        return wave*(1.0+2.0*hmc*wave*(Q/(4.0*np.pi))**2)**2
193
194def GetAsfMean(ElList,Sthl2):
195#   Calculate various scattering factor terms for PDF calcs
196#   ElList: element dictionary contains scattering factor coefficients, etc.
197#   Sthl2: numpy array of sin theta/lambda squared values
198#   returns: mean(f^2), mean(f)^2, mean(compton)
199    sumNoAtoms = 0.0
200    FF = np.zeros_like(Sthl2)
201    FF2 = np.zeros_like(Sthl2)
202    CF = np.zeros_like(Sthl2)
203    for El in ElList:
204        sumNoAtoms += ElList[El]['FormulaNo']
205    for El in ElList:
206        el = ElList[El]
207        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
208        cf = G2elem.ComptonFac(el,Sthl2)
209        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
210        FF2 += ff2*el['FormulaNo']/sumNoAtoms
211        CF += cf*el['FormulaNo']/sumNoAtoms
212    return FF2,FF**2,CF
213
214def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
215    NS = 500
216    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
217    Cgama = np.cos(Gama)[:,np.newaxis]
218    Sgama = np.sin(Gama)[:,np.newaxis]
219    CSgama = 1.0/Sgama
220    Delt = Gama[1]-Gama[0]
221    emc = 7.94e-26
222    Navo = 6.023e23
223    Cth = npcosd(Tth/2.0)
224    CTth = npcosd(Tth)
225    Sth = npcosd(Tth/2.0)
226    STth = npsind(Tth)
227    CSth = 1.0/Sth
228    CSTth = 1.0/STth
229    SCth = 1.0/Cth
230    SCTth = 1.0/CTth
231    if 'Bragg' in Geometry:
232        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
233        Y1 = np.pi
234        Y2 = np.pi/2.0
235        Y3 = 3.*np.pi/8. #3pi/4?
236        W = 1.0/(Sth+np.fabs(Sgama))
237        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
238            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
239        Fac0 = Sth**2*Sgama**2
240        X = Fac0+(Fac0+CTth)**2/2
241        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
242        Z = Cgama**4*Cth**4/2.0
243        E = 2.0*(1.0-a)/(b*Cgama/Cth)
244        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
245        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
246        T1 = np.pi/np.sqrt(F1**2-1.0)
247        T2 = np.pi/np.sqrt(F2**2-1.0)
248        Y4 = T1+T2
249        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
250        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
251        Y7 = (T2-T1)/(F1-F2)
252        YT = F2**2*T2-F1**2*T1
253        Y8 = Y1+YT/(F1-F2)
254        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
255        M = (a**2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E**2*(X*Y7+Y*Y8+Z*Y9))*Cgama
256
257        Q = np.sum(W*M,axis=0)
258        return Q*G*NS/(NS-1.)
259#
260#      cos2th=2.0d*costh^2 - 1.0d
261#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
262#      Y1=3.1415926d
263#      Y2=Y1*0.5d
264#      Y3=Y2*0.75d
265#      for i=1,num_steps-1 do begin
266#         cosgama=double(cos(gama[i]))
267#         singama=double(sin(gama[i]))
268#         cscgama=1.0d / singama
269#
270#         W=1.0d/(sinth+abs(singama))
271#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- \$
272#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
273#
274#         factor0=sinth^2*singama^2
275#         X=factor0+(factor0+cos2th)^2/2.0d
276#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
277#         Z=cosgama^4/2.0d*costh^4
278#         E=2.0d*(1.0-a)/b/cosgama/costh
279#
280#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
281#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
282#
283#         T1=3.14159/sqrt(F1^2-1.0d)
284#         T2=3.14159/sqrt(F2^2-1.0d)
285#         Y4=T1+T2
286#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
287#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
288#         Y7=(T2-T1)/(F1-F2)
289#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
290#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
291#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* \$
292#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
293#
294#         Q=Q+W*M
295#
296#      endfor
297#      Q=double(num_steps)/Double(num_steps-1)*Q*G
298#      end
299    elif 'Cylinder' in Geometry:
300        Q = 0.
301    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
302        EMA = np.exp(-MuT*(1.0-SCTth))
303        Fac1 = (1.-EMA)/(1.0-SCTth)
304        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
305        Fac0 = Cgama/(1-Sgama*SCTth)
306        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
307        Fac0 = Cgama/(1.0+Sgama*SCTth)
308        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
309        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
310        Y = Sgama**3*Cgama*STth*CTth
311        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
312        Fac2 = 1.0/(b*Cgama*STth)
313        U = 2.0*(1.0-a)*Fac2
314        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
315        Mp = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0-Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
316        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
317        Y = -Y
318        Mm = 2.0*np.pi*(a+2.0*(1.0-a)/(2.0+b*(1.0+Sgama)))*(a*X+a*Z/2.0-U*Y+U*(X+Y*V+Z*V**2)/np.sqrt(V**2-1.0)-U*Z*V)
319        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
320        return Q*G*NS/(NS-1.)
321    elif 'Tilting' in Geometry:
322        EMA = np.exp(-MuT*SCth)
323        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
324#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* \$
325#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
326#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* \$
327#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
328        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
329#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
330        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
331#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
332        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
333        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
334        Z = 0.5*Cgama**4*Sth**4
335#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - \$
336#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
337#
338#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
339#
340#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
341#
342        AlP = 2.0+b*(1.0-Cth*Sgama)
343        AlM = 2.0+b*(1.0+Cth*Sgama)
344#          alphaplus = 2.0 + b*(1.0 - costh*singama)
345#          alphaminus = 2.0 + b*(1.0 + costh*singama)
346        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
347        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
348#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
349#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
350        Mp = Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlP)* \
351            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
352            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
353#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + \$
354#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* \$
355#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + \$
356#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - \$
357#                   cosgama*cosgama/b/b*sinth*sinth))
358        Mm =Cgama*(np.pi*a**2*(2.0*X+Y+0.75*Z)+(2.0*np.pi*(1.0-a))*(1.0-a+a*AlM)* \
359            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
360            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
361#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + \$
362#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* \$
363#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + \$
364#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - \$
365#                   cosgama*cosgama/b/b*sinth*sinth))
366        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
367        return Q*G*NS/(NS-1.)
368#       expmutsecth = Exp(-1.0*mut*secth)
369#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
370#       for i=1, num_steps-1 do begin
371#          cosgama=double(cos(gama[i]))
372#          singama=double(sin(gama[i]))
373#          cscgama=1.0d/singama
374#
375#
376#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
377#
378#
379#
380#
381#    ;               print, a,b
382#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
383#          Q=Q+ Wplus*Mplus + Wminus*Mminus
384#      endfor
385#      Q=double(num_steps)/double(num_steps-1)*Q*G
386#   ;   print, min(q), max(q)
387#      end
388
389def MultiScattering(Geometry,ElList,Tth):
390    BN = BD = 0.0
391    Amu = 0.0
392    for El in ElList:
393        el = ElList[El]
394        BN += el['Z']*el['FormulaNo']
395        BD += el['FormulaNo']
396        Amu += el['FormulaNo']*el['mu']
397
398
399def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
400    # returns value(esd) string; nTZ=True for no trailing zeros
401    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
402    #get the 2 significant digits in the esd
403    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
404    #get the number of digits to represent them
405    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
406
407    mdec = lambda esd: -int(math.log10(abs(esd)))
408    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
409    if esd > 0:
410        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
411        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
412        return fmt%(value,int(esd*10**(mdec(esd)+1)))
413    elif esd < 0:
414         return str(round(value,mdec(esd)))
415    else:
416        text = "%F"%(value)
417        if nTZ:
418            return text.rstrip('0')
419        else:
420            return text
421
422
423#GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model
424
425def DoPeakFit(peaks,background,limits,inst,data):
426
427    def backgroundPrint(background,sigback):
428        if background[1]:
429            print 'Background coefficients for',background[0],'function'
430            ptfmt = "%12.5f"
431            ptstr =  'values:'
432            sigstr = 'esds  :'
433            for i,back in enumerate(background[3:]):
434                ptstr += ptfmt % (back)
435                sigstr += ptfmt % (sigback[i+3])
436            print ptstr
437            print sigstr
438        else:
439            print 'Background not refined'
440
441    def instPrint(instVal,siginst,insLabels):
442        print 'Instrument Parameters:'
443        ptfmt = "%12.6f"
444        ptlbls = 'names :'
445        ptstr =  'values:'
446        sigstr = 'esds  :'
447        for i,value in enumerate(instVal):
448            ptlbls += "%s" % (insLabels[i].center(12))
449            ptstr += ptfmt % (value)
450            if siginst[i]:
451                sigstr += ptfmt % (siginst[i])
452            else:
453                sigstr += 12*' '
454        print ptlbls
455        print ptstr
456        print sigstr
457
458    def peaksPrint(peaks,sigpeaks):
459        print 'Peak list='
460
461    begin = time.time()
462    Np = len(peaks[0])
463    DataType = inst[1][0]
464    instVal = inst[1][1:]
465    Insref = inst[2][1:]
466    insLabels = inst[3][1:]
467    Ka2 = False
468    Ioff = 3
469    if len(instVal) == 12:
470        lamratio = instVal[1]/instVal[0]
471        Ka2 = True
472        Ioff = 5
473    insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
474    for peak in peaks:
475        dip = []
476        dip.append(tand(peak[0]/2.0)**2)
477        dip.append(tand(peak[0]/2.0))
478        dip.append(1.0/cosd(peak[0]/2.0))
479        dip.append(tand(peak[0]/2.0))
480        peak.append(dip)
481    B = background[2]
482    bcof = background[3:3+B]
483    Bv = 0
484    if background[1]:
485        Bv = B
486    x,y,w,yc,yb,yd = data               #these are numpy arrays!
487    V = []
488    A = []
489    swobs = 0.0
490    smin = 0.0
491    first = True
492    GoOn = True
493    Go = True
494    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
495        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
496    screenSize = wx.DisplaySize()
497    Size = dlg.GetSize()
498    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
499    try:
500        i = 0
501        for xi in x :
502            Go = dlg.Update(i)[0]
503            if GoOn:
504                GoOn = Go
505            if limits[0] <= xi <= limits[1]:
506                yb[i] = 0.0
507                dp = []
508                for j in range(B):
509                    t = (xi-limits[0])**j
510                    yb[i] += t*bcof[j]
511                    if background[1]:
512                        dp.append(t)
513                yc[i] = yb[i]
514                Iv = 0
515                for j in range(6):
516                    if insref[j]:
517                        dp.append(0.0)
518                        Iv += 1
519                for peak in peaks:
520                    dip = peak[-1]
521                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0)
522                    yc[i] += f[0]*peak[2]
523                    if f[0] > 0.0:
524                        j = 0
525                        if insref[0]:              #U
526                            dp[Bv+j] += f[3]*dip[0]
527                            j += 1
528                        if insref[1]:              #V
529                            dp[Bv+j] += f[3]*dip[1]
530                            j += 1
531                        if insref[2]:              #W
532                            dp[Bv+j] += f[3]
533                            j += 1
534                        if insref[3]:              #X
535                            dp[Bv+j] += f[4]*dip[2]
536                            j += 1
537                        if insref[4]:              #Y
538                            dp[Bv+j] += f[4]*dip[3]
539                            j += 1
540                        if insref[5]:              #SH/L
541                            dp[Bv+j] += f[5]
542                    if Ka2:
543                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
544                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0)
545                       yc[i] += f2[0]*peak[2]*instVal[3]
546                       if f[0] > 0.0:
547                           j = 0
548                           if insref[0]:              #U
549                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
550                               j += 1
551                           if insref[1]:              #V
552                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
553                               j += 1
554                           if insref[2]:              #W
555                               dp[Bv+j] += f2[3]*instVal[3]
556                               j += 1
557                           if insref[3]:              #X
558                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
559                               j += 1
560                           if insref[4]:              #Y
561                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
562                               j += 1
563                           if insref[5]:              #SH/L
564                               dp[Bv+j] += f2[5]*instVal[3]
565                    for j in range(0,Np,2):
566                        if peak[j+1]: dp.append(f[j/2+1])
567                yd[i] = y[i]-yc[i]
568                swobs += w[i]*y[i]**2
569                t2 = w[i]*yd[i]
570                smin += t2*yd[i]
571                if first:
572                    first = False
573                    M = len(dp)
574                    A = np.zeros(shape=(M,M))
575                    V = np.zeros(shape=(M))
576                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
577            i += 1
578    finally:
579        dlg.Destroy()
580    Rwp = smin/swobs
581    Rwp = math.sqrt(Rwp)*100.0
582    norm = np.diag(A)
583    for i,elm in enumerate(norm):
584        if elm <= 0.0:
585            print norm
586            return False,0,0,0,False
587    for i in xrange(len(V)):
588        norm[i] = 1.0/math.sqrt(norm[i])
589        V[i] *= norm[i]
590        a = A[i]
591        for j in xrange(len(V)):
592            a[j] *= norm[i]
593    A = np.transpose(A)
594    for i in xrange(len(V)):
595        a = A[i]
596        for j in xrange(len(V)):
597            a[j] *= norm[i]
598    b = nl.solve(A,V)
599    A = nl.inv(A)
600    sig = np.diag(A)
601    for i in xrange(len(V)):
602        b[i] *= norm[i]
603        sig[i] *= norm[i]*norm[i]
604        sig[i] = math.sqrt(abs(sig[i]))
605    sigback = [0,0,0]
606    for j in range(Bv):
607        background[j+3] += b[j]
608        sigback.append(sig[j])
609    backgroundPrint(background,sigback)
610    k = 0
611    delt = []
612    if Ka2:
613        siginst = [0,0,0,0,0]
614    else:
615        siginst = [0,0,0]
616    for j in range(6):
617        if insref[j]:
618            instVal[j+Ioff] += b[Bv+k]
619            siginst.append(sig[Bv+k])
620            delt.append(b[Bv+k])
621            k += 1
622        else:
623            delt.append(0.0)
624            siginst.append(0.0)
625    delt.append(0.0)                    #dummies for azm
626    siginst.append(0.0)
627    instPrint(instVal,siginst,insLabels)
628    inst[1] = [DataType,]
629    for val in instVal:
630        inst[1].append(val)
631    B = Bv+Iv
632    for peak in peaks:
633        del peak[-1]                        # remove dip from end
634        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
635        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
636        for j in range(0,len(peak[:-1]),2):
637            if peak[j+1]:
638                peak[j] += b[B]
639                B += 1
640        peak[4] += delsig
641        if peak[4] < 0.0:
642            print 'ERROR - negative sigma'
643            return False,0,0,0,False
644        peak[6] += delgam
645        if peak[6] < 0.0:
646            print 'ERROR - negative gamma'
647            return False,0,0,0,False
648    runtime = time.time()-begin
649    data = [x,y,w,yc,yb,yd]
650    return True,smin,Rwp,runtime,GoOn
651
652def CalcPDF(data,inst,xydata):
653    auxPlot = []
654    import copy
655    import scipy.fftpack as ft
656    #subtract backgrounds - if any
657    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
658    if data['Sample Bkg.']['Name']:
659        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
661    if data['Container']['Name']:
663        if data['Container Bkg.']['Name']:
664            xycontainer += (xydata['Container Bkg.'][1][1]+
666        xydata['IofQ'][1][1] += xycontainer
667    #get element data & absorption coeff.
668    ElList = data['ElList']
669    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
670    #Apply angle dependent corrections
671    Tth = xydata['Sample'][1][0]
672    dt = (Tth[1]-Tth[0])
673    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
674    pola = Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
675    auxPlot.append([Tth,pola,'Pola'])
676    xydata['IofQ'][1][1] /= pola
677    xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
678    XY = xydata['IofQ'][1]
679    #convert to Q
680    hc = 12.397639
681    if 'Lam' in inst:
682        wave = inst['Lam']
683    else:
684        wave = inst['Lam1']
685    keV = hc/wave
686    minQ = npT2q(Tth[0],wave)
687    maxQ = npT2q(Tth[-1],wave)
688    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
689    dq = Qpoints[1]-Qpoints[0]
690    XY[0] = npT2q(XY[0],wave)
691    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear'))
692
693    qLimits = data['QScaleLim']
694    minQ = np.searchsorted(Qpoints,qLimits[0])
695    maxQ = np.searchsorted(Qpoints,qLimits[1])
696    newdata = []
697    xydata['IofQ'][1][0] = Qpoints
698    xydata['IofQ'][1][1] = Qdata
699    for item in xydata['IofQ'][1]:
700        newdata.append(item[:maxQ])
701    xydata['IofQ'][1] = newdata
702
703
704    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
705    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
706    Q = xydata['SofQ'][1][0]
707    ruland = Ruland(data['Ruland'],wave,Q,CF)
708    auxPlot.append([Q,ruland,'Ruland'])
709    CF *= ruland
710    auxPlot.append([Q,CF,'Compton*Ruland'])
711    scale =       np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
712    xydata['SofQ'][1][1] *= scale
713    xydata['SofQ'][1][1] -= CF
714    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
715
716    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
717    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
718    xydata['FofQ'][1][1] *= np.sin(np.pi*(qLimits[1]-Q)/(2.0*qLimits[1]))
719
720    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
721    nR = len(xydata['GofR'][1][1])
722    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
723    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
724
725
726    return auxPlot
727
Note: See TracBrowser for help on using the repository browser.