1 | #GSASII powder calculation module |
---|
2 | ########### SVN repository information ################### |
---|
3 | # $Date: 2011-04-20 13:09:53 -0500 (Wed, 20 Apr 2011) $ |
---|
4 | # $Author: vondreele $ |
---|
5 | # $Revision: 267 $ |
---|
6 | # $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIpwd.py $ |
---|
7 | # $Id: GSASIIpwd.py 267 2011-04-20 18:09:53Z vondreele $ |
---|
8 | ########### SVN repository information ################### |
---|
9 | import sys |
---|
10 | import math |
---|
11 | import wx |
---|
12 | import time |
---|
13 | import numpy as np |
---|
14 | import scipy as sp |
---|
15 | import numpy.linalg as nl |
---|
16 | import scipy.interpolate as si |
---|
17 | import GSASIIpath |
---|
18 | import pypowder as pyp #assumes path has been amended to include correctr bin directory |
---|
19 | import GSASIIplot as G2plt |
---|
20 | import GSASIIlattice as G2lat |
---|
21 | import GSASIIElem as G2elem |
---|
22 | import GSASIIgrid as G2gd |
---|
23 | |
---|
24 | # trig functions in degrees |
---|
25 | sind = lambda x: math.sin(x*math.pi/180.) |
---|
26 | asind = lambda x: 180.*math.asin(x)/math.pi |
---|
27 | tand = lambda x: math.tan(x*math.pi/180.) |
---|
28 | atand = lambda x: 180.*math.atan(x)/math.pi |
---|
29 | atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi |
---|
30 | cosd = lambda x: math.cos(x*math.pi/180.) |
---|
31 | acosd = lambda x: 180.*math.acos(x)/math.pi |
---|
32 | rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p) |
---|
33 | #numpy versions |
---|
34 | npsind = lambda x: np.sin(x*np.pi/180.) |
---|
35 | npasind = lambda x: 180.*np.arcsin(x)/math.pi |
---|
36 | npcosd = lambda x: np.cos(x*math.pi/180.) |
---|
37 | nptand = lambda x: np.tan(x*math.pi/180.) |
---|
38 | npatand = lambda x: 180.*np.arctan(x)/np.pi |
---|
39 | npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
40 | npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave |
---|
41 | npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave) |
---|
42 | |
---|
43 | def 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 | |
---|
68 | def 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 | |
---|
85 | def 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 | |
---|
115 | def 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 | |
---|
160 | def 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 ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ |
---|
169 | (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 |
---|
170 | |
---|
171 | def Oblique(ObCoeff,Tth): |
---|
172 | if ObCoeff: |
---|
173 | return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth))) |
---|
174 | else: |
---|
175 | return 1.0 |
---|
176 | |
---|
177 | def Ruland(RulCoff,wave,Q,Compton): |
---|
178 | C = 2.9978e8 |
---|
179 | D = 1.5e-3 |
---|
180 | hmc = 0.024262734687 |
---|
181 | sinth2 = (Q*wave/(4.0*np.pi))**2 |
---|
182 | dlam = (wave**2)*Compton*Q/C |
---|
183 | dlam_c = 2.0*hmc*sinth2-D*wave**2 |
---|
184 | return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2)) |
---|
185 | |
---|
186 | def LorchWeight(Q): |
---|
187 | return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1])) |
---|
188 | |
---|
189 | def GetAsfMean(ElList,Sthl2): |
---|
190 | # Calculate various scattering factor terms for PDF calcs |
---|
191 | # ElList: element dictionary contains scattering factor coefficients, etc. |
---|
192 | # Sthl2: numpy array of sin theta/lambda squared values |
---|
193 | # returns: mean(f^2), mean(f)^2, mean(compton) |
---|
194 | sumNoAtoms = 0.0 |
---|
195 | FF = np.zeros_like(Sthl2) |
---|
196 | FF2 = np.zeros_like(Sthl2) |
---|
197 | CF = np.zeros_like(Sthl2) |
---|
198 | for El in ElList: |
---|
199 | sumNoAtoms += ElList[El]['FormulaNo'] |
---|
200 | for El in ElList: |
---|
201 | el = ElList[El] |
---|
202 | ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2 |
---|
203 | cf = G2elem.ComptonFac(el,Sthl2) |
---|
204 | FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms |
---|
205 | FF2 += ff2*el['FormulaNo']/sumNoAtoms |
---|
206 | CF += cf*el['FormulaNo']/sumNoAtoms |
---|
207 | return FF2,FF**2,CF |
---|
208 | |
---|
209 | def GetNumDensity(ElList,Vol): |
---|
210 | sumNoAtoms = 0.0 |
---|
211 | for El in ElList: |
---|
212 | sumNoAtoms += ElList[El]['FormulaNo'] |
---|
213 | return sumNoAtoms/Vol |
---|
214 | |
---|
215 | def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01): |
---|
216 | NS = 500 |
---|
217 | Gama = np.linspace(0.,np.pi/2.,NS,False)[1:] |
---|
218 | Cgama = np.cos(Gama)[:,np.newaxis] |
---|
219 | Sgama = np.sin(Gama)[:,np.newaxis] |
---|
220 | CSgama = 1.0/Sgama |
---|
221 | Delt = Gama[1]-Gama[0] |
---|
222 | emc = 7.94e-26 |
---|
223 | Navo = 6.023e23 |
---|
224 | Cth = npcosd(Tth/2.0) |
---|
225 | CTth = npcosd(Tth) |
---|
226 | Sth = npcosd(Tth/2.0) |
---|
227 | STth = npsind(Tth) |
---|
228 | CSth = 1.0/Sth |
---|
229 | CSTth = 1.0/STth |
---|
230 | SCth = 1.0/Cth |
---|
231 | SCTth = 1.0/CTth |
---|
232 | if 'Bragg' in Geometry: |
---|
233 | G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth))) |
---|
234 | Y1 = np.pi |
---|
235 | Y2 = np.pi/2.0 |
---|
236 | Y3 = 3.*np.pi/8. #3pi/4? |
---|
237 | W = 1.0/(Sth+np.fabs(Sgama)) |
---|
238 | W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))- |
---|
239 | (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2) |
---|
240 | Fac0 = Sth**2*Sgama**2 |
---|
241 | X = Fac0+(Fac0+CTth)**2/2 |
---|
242 | Y = Cgama**2*Cth**2*(1.0-Fac0-CTth) |
---|
243 | Z = Cgama**4*Cth**4/2.0 |
---|
244 | E = 2.0*(1.0-a)/(b*Cgama/Cth) |
---|
245 | F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1 |
---|
246 | F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama) |
---|
247 | T1 = np.pi/np.sqrt(F1**2-1.0) |
---|
248 | T2 = np.pi/np.sqrt(F2**2-1.0) |
---|
249 | Y4 = T1+T2 |
---|
250 | Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2) |
---|
251 | Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3) |
---|
252 | Y7 = (T2-T1)/(F1-F2) |
---|
253 | YT = F2**2*T2-F1**2*T1 |
---|
254 | Y8 = Y1+YT/(F1-F2) |
---|
255 | Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2) |
---|
256 | 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 |
---|
257 | |
---|
258 | Q = np.sum(W*M,axis=0) |
---|
259 | return Q*G*NS/(NS-1.) |
---|
260 | # |
---|
261 | # cos2th=2.0d*costh^2 - 1.0d |
---|
262 | # G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth)) |
---|
263 | # Y1=3.1415926d |
---|
264 | # Y2=Y1*0.5d |
---|
265 | # Y3=Y2*0.75d |
---|
266 | # for i=1,num_steps-1 do begin |
---|
267 | # cosgama=double(cos(gama[i])) |
---|
268 | # singama=double(sin(gama[i])) |
---|
269 | # cscgama=1.0d / singama |
---|
270 | # |
---|
271 | # W=1.0d/(sinth+abs(singama)) |
---|
272 | # W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $ |
---|
273 | # (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2) |
---|
274 | # |
---|
275 | # factor0=sinth^2*singama^2 |
---|
276 | # X=factor0+(factor0+cos2th)^2/2.0d |
---|
277 | # Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2 |
---|
278 | # Z=cosgama^4/2.0d*costh^4 |
---|
279 | # E=2.0d*(1.0-a)/b/cosgama/costh |
---|
280 | # |
---|
281 | # F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh |
---|
282 | # F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh |
---|
283 | # |
---|
284 | # T1=3.14159/sqrt(F1^2-1.0d) |
---|
285 | # T2=3.14159/sqrt(F2^2-1.0d) |
---|
286 | # Y4=T1+T2 |
---|
287 | # Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2) |
---|
288 | # Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3) |
---|
289 | # Y7=(T2-T1)/(F1-F2) |
---|
290 | # Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2) |
---|
291 | # Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2) |
---|
292 | # M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $ |
---|
293 | # (X*Y7+Y*Y8+Z*Y9))*cosgama |
---|
294 | # |
---|
295 | # Q=Q+W*M |
---|
296 | # |
---|
297 | # endfor |
---|
298 | # Q=double(num_steps)/Double(num_steps-1)*Q*G |
---|
299 | # end |
---|
300 | elif 'Cylinder' in Geometry: |
---|
301 | Q = 0. |
---|
302 | elif 'Fixed' in Geometry: #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections |
---|
303 | EMA = np.exp(-MuT*(1.0-SCTth)) |
---|
304 | Fac1 = (1.-EMA)/(1.0-SCTth) |
---|
305 | G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1) |
---|
306 | Fac0 = Cgama/(1-Sgama*SCTth) |
---|
307 | Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0)) |
---|
308 | Fac0 = Cgama/(1.0+Sgama*SCTth) |
---|
309 | Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0)) |
---|
310 | X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0 |
---|
311 | Y = Sgama**3*Cgama*STth*CTth |
---|
312 | Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0 |
---|
313 | Fac2 = 1.0/(b*Cgama*STth) |
---|
314 | U = 2.0*(1.0-a)*Fac2 |
---|
315 | V = (2.0+b*(1.0-CTth*Sgama))*Fac2 |
---|
316 | 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) |
---|
317 | V = (2.0+b*(1.0+CTth*Sgama))*Fac2 |
---|
318 | Y = -Y |
---|
319 | 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) |
---|
320 | Q = np.sum(Wp*Mp+Wm*Mm,axis=0) |
---|
321 | return Q*G*NS/(NS-1.) |
---|
322 | elif 'Tilting' in Geometry: |
---|
323 | EMA = np.exp(-MuT*SCth) |
---|
324 | G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA) |
---|
325 | # Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $ |
---|
326 | # (Exp(-1.0*mut*cscgama) - expmutsecth) |
---|
327 | # Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $ |
---|
328 | # expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama)) |
---|
329 | Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA) |
---|
330 | # Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA) |
---|
331 | Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama)) |
---|
332 | # Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama)) |
---|
333 | X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0) |
---|
334 | Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2 |
---|
335 | Z = 0.5*Cgama**4*Sth**4 |
---|
336 | # X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $ |
---|
337 | # 4.0*sinth*sinth*cosgama*cosgama) +1.0d) |
---|
338 | # |
---|
339 | # Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth |
---|
340 | # |
---|
341 | # Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4) |
---|
342 | # |
---|
343 | AlP = 2.0+b*(1.0-Cth*Sgama) |
---|
344 | AlM = 2.0+b*(1.0+Cth*Sgama) |
---|
345 | # alphaplus = 2.0 + b*(1.0 - costh*singama) |
---|
346 | # alphaminus = 2.0 + b*(1.0 + costh*singama) |
---|
347 | BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2)) |
---|
348 | BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2)) |
---|
349 | # betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth)) |
---|
350 | # betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth)) |
---|
351 | 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)* \ |
---|
352 | (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+ |
---|
353 | 2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2)) |
---|
354 | # Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $ |
---|
355 | # (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $ |
---|
356 | # (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $ |
---|
357 | # 2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $ |
---|
358 | # cosgama*cosgama/b/b*sinth*sinth)) |
---|
359 | 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)* \ |
---|
360 | (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+ |
---|
361 | 2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2)) |
---|
362 | # Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $ |
---|
363 | # (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $ |
---|
364 | # (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $ |
---|
365 | # 2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $ |
---|
366 | # cosgama*cosgama/b/b*sinth*sinth)) |
---|
367 | Q = np.sum(Wp*Mp+Wm*Mm,axis=0) |
---|
368 | return Q*G*NS/(NS-1.) |
---|
369 | # expmutsecth = Exp(-1.0*mut*secth) |
---|
370 | # G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth |
---|
371 | # for i=1, num_steps-1 do begin |
---|
372 | # cosgama=double(cos(gama[i])) |
---|
373 | # singama=double(sin(gama[i])) |
---|
374 | # cscgama=1.0d/singama |
---|
375 | # |
---|
376 | # |
---|
377 | # ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus) |
---|
378 | # |
---|
379 | # |
---|
380 | # |
---|
381 | # |
---|
382 | # ; print, a,b |
---|
383 | # ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus) |
---|
384 | # Q=Q+ Wplus*Mplus + Wminus*Mminus |
---|
385 | # endfor |
---|
386 | # Q=double(num_steps)/double(num_steps-1)*Q*G |
---|
387 | # ; print, min(q), max(q) |
---|
388 | # end |
---|
389 | |
---|
390 | def MultiScattering(Geometry,ElList,Tth): |
---|
391 | BN = BD = 0.0 |
---|
392 | Amu = 0.0 |
---|
393 | for El in ElList: |
---|
394 | el = ElList[El] |
---|
395 | BN += el['Z']*el['FormulaNo'] |
---|
396 | BD += el['FormulaNo'] |
---|
397 | Amu += el['FormulaNo']*el['mu'] |
---|
398 | |
---|
399 | |
---|
400 | def ValEsd(value,esd=0,nTZ=False): #NOT complete - don't use |
---|
401 | # returns value(esd) string; nTZ=True for no trailing zeros |
---|
402 | # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal |
---|
403 | #get the 2 significant digits in the esd |
---|
404 | edig = lambda esd: int(round(10**(math.log10(esd) % 1+1))) |
---|
405 | #get the number of digits to represent them |
---|
406 | epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd))) |
---|
407 | |
---|
408 | mdec = lambda esd: -int(math.log10(abs(esd))) |
---|
409 | ndec = lambda esd: int(1.545-math.log10(abs(esd))) |
---|
410 | if esd > 0: |
---|
411 | fmt = '"%.'+str(ndec(esd))+'f(%d)"' |
---|
412 | print fmt,ndec(esd),esd*10**(mdec(esd)+1) |
---|
413 | return fmt%(value,int(esd*10**(mdec(esd)+1))) |
---|
414 | elif esd < 0: |
---|
415 | return str(round(value,mdec(esd))) |
---|
416 | else: |
---|
417 | text = "%F"%(value) |
---|
418 | if nTZ: |
---|
419 | return text.rstrip('0') |
---|
420 | else: |
---|
421 | return text |
---|
422 | |
---|
423 | |
---|
424 | #GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model |
---|
425 | |
---|
426 | def DoPeakFit(peaks,background,limits,inst,data): |
---|
427 | |
---|
428 | def backgroundPrint(background,sigback): |
---|
429 | if background[1]: |
---|
430 | print 'Background coefficients for',background[0],'function' |
---|
431 | ptfmt = "%12.5f" |
---|
432 | ptstr = 'values:' |
---|
433 | sigstr = 'esds :' |
---|
434 | for i,back in enumerate(background[3:]): |
---|
435 | ptstr += ptfmt % (back) |
---|
436 | sigstr += ptfmt % (sigback[i+3]) |
---|
437 | print ptstr |
---|
438 | print sigstr |
---|
439 | else: |
---|
440 | print 'Background not refined' |
---|
441 | |
---|
442 | def instPrint(instVal,siginst,insLabels): |
---|
443 | print 'Instrument Parameters:' |
---|
444 | ptfmt = "%12.6f" |
---|
445 | ptlbls = 'names :' |
---|
446 | ptstr = 'values:' |
---|
447 | sigstr = 'esds :' |
---|
448 | for i,value in enumerate(instVal): |
---|
449 | ptlbls += "%s" % (insLabels[i].center(12)) |
---|
450 | ptstr += ptfmt % (value) |
---|
451 | if siginst[i]: |
---|
452 | sigstr += ptfmt % (siginst[i]) |
---|
453 | else: |
---|
454 | sigstr += 12*' ' |
---|
455 | print ptlbls |
---|
456 | print ptstr |
---|
457 | print sigstr |
---|
458 | |
---|
459 | def peaksPrint(peaks,sigpeaks): |
---|
460 | print 'Peak list=' |
---|
461 | |
---|
462 | begin = time.time() |
---|
463 | Np = len(peaks[0]) |
---|
464 | DataType = inst[1][0] |
---|
465 | instVal = inst[1][1:] |
---|
466 | Insref = inst[2][1:] |
---|
467 | insLabels = inst[3][1:] |
---|
468 | Ka2 = False |
---|
469 | Ioff = 3 |
---|
470 | if len(instVal) == 12: |
---|
471 | lamratio = instVal[1]/instVal[0] |
---|
472 | Ka2 = True |
---|
473 | Ioff = 5 |
---|
474 | insref = Insref[len(Insref)-7:-1] #just U,V,W,X,Y,SH/L |
---|
475 | for peak in peaks: |
---|
476 | dip = [] |
---|
477 | dip.append(tand(peak[0]/2.0)**2) |
---|
478 | dip.append(tand(peak[0]/2.0)) |
---|
479 | dip.append(1.0/cosd(peak[0]/2.0)) |
---|
480 | dip.append(tand(peak[0]/2.0)) |
---|
481 | peak.append(dip) |
---|
482 | B = background[2] |
---|
483 | bcof = background[3:3+B] |
---|
484 | Bv = 0 |
---|
485 | if background[1]: |
---|
486 | Bv = B |
---|
487 | x,y,w,yc,yb,yd = data #these are numpy arrays! |
---|
488 | V = [] |
---|
489 | A = [] |
---|
490 | swobs = 0.0 |
---|
491 | smin = 0.0 |
---|
492 | first = True |
---|
493 | GoOn = True |
---|
494 | Go = True |
---|
495 | dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \ |
---|
496 | style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT) |
---|
497 | screenSize = wx.DisplaySize() |
---|
498 | Size = dlg.GetSize() |
---|
499 | dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0)) |
---|
500 | try: |
---|
501 | i = 0 |
---|
502 | for xi in x : |
---|
503 | Go = dlg.Update(i)[0] |
---|
504 | if GoOn: |
---|
505 | GoOn = Go |
---|
506 | if limits[0] <= xi <= limits[1]: |
---|
507 | yb[i] = 0.0 |
---|
508 | dp = [] |
---|
509 | for j in range(B): |
---|
510 | t = (xi-limits[0])**j |
---|
511 | yb[i] += t*bcof[j] |
---|
512 | if background[1]: |
---|
513 | dp.append(t) |
---|
514 | yc[i] = yb[i] |
---|
515 | Iv = 0 |
---|
516 | for j in range(6): |
---|
517 | if insref[j]: |
---|
518 | dp.append(0.0) |
---|
519 | Iv += 1 |
---|
520 | for peak in peaks: |
---|
521 | dip = peak[-1] |
---|
522 | f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],instVal[-2],0.0) |
---|
523 | yc[i] += f[0]*peak[2] |
---|
524 | if f[0] > 0.0: |
---|
525 | j = 0 |
---|
526 | if insref[0]: #U |
---|
527 | dp[Bv+j] += f[3]*dip[0] |
---|
528 | j += 1 |
---|
529 | if insref[1]: #V |
---|
530 | dp[Bv+j] += f[3]*dip[1] |
---|
531 | j += 1 |
---|
532 | if insref[2]: #W |
---|
533 | dp[Bv+j] += f[3] |
---|
534 | j += 1 |
---|
535 | if insref[3]: #X |
---|
536 | dp[Bv+j] += f[4]*dip[2] |
---|
537 | j += 1 |
---|
538 | if insref[4]: #Y |
---|
539 | dp[Bv+j] += f[4]*dip[3] |
---|
540 | j += 1 |
---|
541 | if insref[5]: #SH/L |
---|
542 | dp[Bv+j] += f[5] |
---|
543 | if Ka2: |
---|
544 | pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0)) |
---|
545 | f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],instVal[-2],0.0) |
---|
546 | yc[i] += f2[0]*peak[2]*instVal[3] |
---|
547 | if f[0] > 0.0: |
---|
548 | j = 0 |
---|
549 | if insref[0]: #U |
---|
550 | dp[Bv+j] += f2[3]*dip[0]*instVal[3] |
---|
551 | j += 1 |
---|
552 | if insref[1]: #V |
---|
553 | dp[Bv+j] += f2[3]*dip[1]*instVal[3] |
---|
554 | j += 1 |
---|
555 | if insref[2]: #W |
---|
556 | dp[Bv+j] += f2[3]*instVal[3] |
---|
557 | j += 1 |
---|
558 | if insref[3]: #X |
---|
559 | dp[Bv+j] += f2[4]*dip[2]*instVal[3] |
---|
560 | j += 1 |
---|
561 | if insref[4]: #Y |
---|
562 | dp[Bv+j] += f2[4]*dip[3]*instVal[3] |
---|
563 | j += 1 |
---|
564 | if insref[5]: #SH/L |
---|
565 | dp[Bv+j] += f2[5]*instVal[3] |
---|
566 | for j in range(0,Np,2): |
---|
567 | if peak[j+1]: dp.append(f[j/2+1]) |
---|
568 | yd[i] = y[i]-yc[i] |
---|
569 | swobs += w[i]*y[i]**2 |
---|
570 | t2 = w[i]*yd[i] |
---|
571 | smin += t2*yd[i] |
---|
572 | if first: |
---|
573 | first = False |
---|
574 | M = len(dp) |
---|
575 | A = np.zeros(shape=(M,M)) |
---|
576 | V = np.zeros(shape=(M)) |
---|
577 | A,V = pyp.buildmv(t2,w[i],M,dp,A,V) |
---|
578 | i += 1 |
---|
579 | finally: |
---|
580 | dlg.Destroy() |
---|
581 | Rwp = smin/swobs |
---|
582 | Rwp = math.sqrt(Rwp)*100.0 |
---|
583 | norm = np.diag(A) |
---|
584 | for i,elm in enumerate(norm): |
---|
585 | if elm <= 0.0: |
---|
586 | print norm |
---|
587 | return False,0,0,0,False |
---|
588 | for i in xrange(len(V)): |
---|
589 | norm[i] = 1.0/math.sqrt(norm[i]) |
---|
590 | V[i] *= norm[i] |
---|
591 | a = A[i] |
---|
592 | for j in xrange(len(V)): |
---|
593 | a[j] *= norm[i] |
---|
594 | A = np.transpose(A) |
---|
595 | for i in xrange(len(V)): |
---|
596 | a = A[i] |
---|
597 | for j in xrange(len(V)): |
---|
598 | a[j] *= norm[i] |
---|
599 | b = nl.solve(A,V) |
---|
600 | A = nl.inv(A) |
---|
601 | sig = np.diag(A) |
---|
602 | for i in xrange(len(V)): |
---|
603 | b[i] *= norm[i] |
---|
604 | sig[i] *= norm[i]*norm[i] |
---|
605 | sig[i] = math.sqrt(abs(sig[i])) |
---|
606 | sigback = [0,0,0] |
---|
607 | for j in range(Bv): |
---|
608 | background[j+3] += b[j] |
---|
609 | sigback.append(sig[j]) |
---|
610 | backgroundPrint(background,sigback) |
---|
611 | k = 0 |
---|
612 | delt = [] |
---|
613 | if Ka2: |
---|
614 | siginst = [0,0,0,0,0] |
---|
615 | else: |
---|
616 | siginst = [0,0,0] |
---|
617 | for j in range(6): |
---|
618 | if insref[j]: |
---|
619 | instVal[j+Ioff] += b[Bv+k] |
---|
620 | siginst.append(sig[Bv+k]) |
---|
621 | delt.append(b[Bv+k]) |
---|
622 | k += 1 |
---|
623 | else: |
---|
624 | delt.append(0.0) |
---|
625 | siginst.append(0.0) |
---|
626 | delt.append(0.0) #dummies for azm |
---|
627 | siginst.append(0.0) |
---|
628 | instPrint(instVal,siginst,insLabels) |
---|
629 | inst[1] = [DataType,] |
---|
630 | for val in instVal: |
---|
631 | inst[1].append(val) |
---|
632 | B = Bv+Iv |
---|
633 | for peak in peaks: |
---|
634 | del peak[-1] # remove dip from end |
---|
635 | delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2] |
---|
636 | delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0) |
---|
637 | for j in range(0,len(peak[:-1]),2): |
---|
638 | if peak[j+1]: |
---|
639 | peak[j] += b[B] |
---|
640 | B += 1 |
---|
641 | peak[4] += delsig |
---|
642 | if peak[4] < 0.0: |
---|
643 | print 'ERROR - negative sigma' |
---|
644 | return False,0,0,0,False |
---|
645 | peak[6] += delgam |
---|
646 | if peak[6] < 0.0: |
---|
647 | print 'ERROR - negative gamma' |
---|
648 | return False,0,0,0,False |
---|
649 | runtime = time.time()-begin |
---|
650 | data = [x,y,w,yc,yb,yd] |
---|
651 | return True,smin,Rwp,runtime,GoOn |
---|
652 | |
---|
653 | def CalcPDF(data,inst,xydata): |
---|
654 | auxPlot = [] |
---|
655 | import copy |
---|
656 | import scipy.fftpack as ft |
---|
657 | #subtract backgrounds - if any |
---|
658 | xydata['IofQ'] = copy.deepcopy(xydata['Sample']) |
---|
659 | if data['Sample Bkg.']['Name']: |
---|
660 | xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+ |
---|
661 | data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult'] |
---|
662 | if data['Container']['Name']: |
---|
663 | xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult'] |
---|
664 | if data['Container Bkg.']['Name']: |
---|
665 | xycontainer += (xydata['Container Bkg.'][1][1]+ |
---|
666 | data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult'] |
---|
667 | xydata['IofQ'][1][1] += xycontainer |
---|
668 | #get element data & absorption coeff. |
---|
669 | ElList = data['ElList'] |
---|
670 | Abs = G2lat.CellAbsorption(ElList,data['Form Vol']) |
---|
671 | #Apply angle dependent corrections |
---|
672 | Tth = xydata['Sample'][1][0] |
---|
673 | dt = (Tth[1]-Tth[0]) |
---|
674 | xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth) |
---|
675 | xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth']) |
---|
676 | if data['DetType'] == 'Image plate': |
---|
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')) #only OK for scipy 0.9! |
---|
692 | T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0) #OK for scipy 0.8 |
---|
693 | Qdata = T(Qpoints) |
---|
694 | |
---|
695 | qLimits = data['QScaleLim'] |
---|
696 | minQ = np.searchsorted(Qpoints,qLimits[0]) |
---|
697 | maxQ = np.searchsorted(Qpoints,qLimits[1]) |
---|
698 | newdata = [] |
---|
699 | xydata['IofQ'][1][0] = Qpoints |
---|
700 | xydata['IofQ'][1][1] = Qdata |
---|
701 | for item in xydata['IofQ'][1]: |
---|
702 | newdata.append(item[:maxQ]) |
---|
703 | xydata['IofQ'][1] = newdata |
---|
704 | |
---|
705 | |
---|
706 | xydata['SofQ'] = copy.deepcopy(xydata['IofQ']) |
---|
707 | FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2) #these are <f^2>,<f>^2,Cf |
---|
708 | Q = xydata['SofQ'][1][0] |
---|
709 | ruland = Ruland(data['Ruland'],wave,Q,CF) |
---|
710 | # auxPlot.append([Q,ruland,'Ruland']) |
---|
711 | CF *= ruland |
---|
712 | # auxPlot.append([Q,CF,'CF-Corr']) |
---|
713 | scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) |
---|
714 | xydata['SofQ'][1][1] *= scale |
---|
715 | xydata['SofQ'][1][1] -= CF |
---|
716 | xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF |
---|
717 | scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) |
---|
718 | xydata['SofQ'][1][1] *= scale |
---|
719 | |
---|
720 | xydata['FofQ'] = copy.deepcopy(xydata['SofQ']) |
---|
721 | xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0) |
---|
722 | if data['Lorch']: |
---|
723 | xydata['FofQ'][1][1] *= LorchWeight(Q) |
---|
724 | |
---|
725 | xydata['GofR'] = copy.deepcopy(xydata['FofQ']) |
---|
726 | nR = len(xydata['GofR'][1][1]) |
---|
727 | xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR]) |
---|
728 | xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1] |
---|
729 | |
---|
730 | |
---|
731 | return auxPlot |
---|
732 | |
---|