source: trunk/GSASIIpwd.py @ 275

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

use interp1d instead of griddata
File size: 28.3 KB
Line 
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 ###################
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 ((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   
171def 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               
177def 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   
186def LorchWeight(Q):
187    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
188           
189def 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   
209def GetNumDensity(ElList,Vol):
210    sumNoAtoms = 0.0
211    for El in ElList:
212        sumNoAtoms += ElList[El]['FormulaNo']
213    return sumNoAtoms/Vol
214           
215def 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
390def 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
400def 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
426def 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
653def 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       
Note: See TracBrowser for help on using the repository browser.