source: trunk/GSASIIpeak.py @ 267

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

removed azm rotation - didn't work
fix azimuthal polarization

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