source: trunk/GSASIIpwd.py @ 271

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

make g77 versions of fortran for binwin2.7 - faster than gfortran!
fix to histogram2d.for
GSASII.py - allow 2D offset in multiplots, fix for large file lists
GSASIIimgGUI.py - more stuff in save/load controls
GSASIIplot.py - 2D offset in multiplots, fix bad behavior in structure drawings
GSASIIpwd.py - more mods to PDF calcs

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