source: trunk/GSASIIpwd.py @ 299

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

Add progress bar to peak fit
File size: 36.3 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2011-04-20 13:09:53 -0500 (Wed, 20 Apr 2011) $
6# $Author: vondreele $
7# $Revision: 267 $
8# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 267 2011-04-20 18:09:53Z vondreele $
10########### SVN repository information ###################
11import sys
12import math
13import wx
14import time
15
16import numpy as np
17import scipy as sp
18import numpy.linalg as nl
19from numpy.fft import ifft, fft, fftshift
20import scipy.interpolate as si
21import scipy.stats as st
22import scipy.optimize as so
23
24import GSASIIpath
25#import pypowder as pyp
26import GSASIIplot as G2plt
27import GSASIIlattice as G2lat
28import GSASIIElem as G2elem
29import GSASIIgrid as G2gd
30import GSASIIIO as G2IO
31
32# trig functions in degrees
33sind = lambda x: math.sin(x*math.pi/180.)
34asind = lambda x: 180.*math.asin(x)/math.pi
35tand = lambda x: math.tan(x*math.pi/180.)
36atand = lambda x: 180.*math.atan(x)/math.pi
37atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
38cosd = lambda x: math.cos(x*math.pi/180.)
39acosd = lambda x: 180.*math.acos(x)/math.pi
40rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
41#numpy versions
42npsind = lambda x: np.sin(x*np.pi/180.)
43npasind = lambda x: 180.*np.arcsin(x)/math.pi
44npcosd = lambda x: np.cos(x*math.pi/180.)
45npacosd = lambda x: 180.*np.arccos(x)/math.pi
46nptand = lambda x: np.tan(x*math.pi/180.)
47npatand = lambda x: 180.*np.arctan(x)/np.pi
48npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
49npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
50npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
51   
52       
53def factorize(num):
54    ''' Provide prime number factors for integer num
55    Returns dictionary of prime factors (keys) & power for each (data)
56    '''
57    factors = {}
58    orig = num
59
60    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
61    i, sqi = 2, 4
62    while sqi <= num:
63        while not num%i:
64            num /= i
65            factors[i] = factors.get(i, 0) + 1
66
67        sqi += 2*i + 1
68        i += 1
69
70    if num != 1 and num != orig:
71        factors[num] = factors.get(num, 0) + 1
72
73    if factors:
74        return factors
75    else:
76        return {num:1}          #a prime number!
77           
78def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
79    ''' Provide list of optimal data sizes for FFT calculations
80    Input:
81        nmin: minimum data size >= 1
82        nmax: maximum data size > nmin
83        thresh: maximum prime factor allowed
84    Returns:
85        list of data sizes where the maximum prime factor is < thresh
86    ''' 
87    plist = []
88    nmin = max(1,nmin)
89    nmax = max(nmin+1,nmax)
90    for p in range(nmin,nmax):
91        if max(factorize(p).keys()) < thresh:
92            plist.append(p)
93    return plist
94
95def ValuesOut(parmdict, varylist):
96    '''Use before call to leastsq to setup list of values for the parameters
97    in parmdict, as selected by key in varylist'''
98    return [parmdict[key] for key in varylist] 
99   
100def ValuesIn(parmdict, varylist, values):
101    ''' Use after call to leastsq to update the parameter dictionary with
102    values corresponding to keys in varylist'''
103    parmdict.update(zip(varylist,values))
104   
105def Transmission(Geometry,Abs,Diam):
106#Calculate sample transmission
107#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
108#   Abs: absorption coeff in cm-1
109#   Diam: sample thickness/diameter in mm
110    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
111        MuR = Abs*Diam/20.0
112        if MuR <= 3.0:
113            T0 = 16/(3.*math.pi)
114            T1 = -0.045780
115            T2 = -0.02489
116            T3 = 0.003045
117            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
118            if T < -20.:
119                return 2.06e-9
120            else:
121                return math.exp(T)
122        else:
123            T1 = 1.433902
124            T2 = 0.013869+0.337894
125            T3 = 1.933433+1.163198
126            T4 = 0.044365-0.04259
127            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
128            return T/100.
129    elif 'plate' in Geometry:
130        MuR = Abs*Diam/10.
131        return math.exp(-MuR)
132    elif 'Bragg' in Geometry:
133        return 0.0
134
135def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
136#Calculate sample absorption
137#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
138#   Abs: absorption coeff in cm-1
139#   Diam: sample thickness/diameter in mm
140#   Tth: 2-theta scattering angle - can be numpy array
141#   Phi: flat plate tilt angle - future
142#   Psi: flat plate tilt axis - future
143    Sth2 = npsind(Tth/2.0)**2
144    Cth2 = 1.-Sth2
145    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
146        MuR = Abs*Diam/20.0
147        if MuR < 3.0:
148            T0 = 16.0/(3*np.pi)
149            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
150                0.109561*np.sqrt(Sth2)-26.04556
151            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
152                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
153            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
154            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
155            return np.exp(Trns)
156        else:
157            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
158                10.02088*Sth2**3-3.36778*Sth2**4
159            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
160                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
161            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
162                0.13576*np.sqrt(Sth2)+1.163198
163            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
164            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
165            return Trns/100.
166    elif 'Bragg' in Geometry:
167        return 1.0
168    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
169        # and only defined for 2theta < 90
170        MuR = Abs*Diam/10.0
171        T1 = np.exp(-MuR)
172        T2 = np.exp(-MuR/npcosd(Tth))
173        Tb = MuR-MuR/npcosd(Tth)
174        return (T2-T1)/Tb
175    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
176        MuR = Abs*Diam/10.0
177        cth = npcosd(Tth/2.0)
178        return np.exp(-MuR/cth)/cth
179       
180def Polarization(Pola,Tth,Azm=0.0):
181#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
182#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
183#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
184#   Tth: 2-theta scattering angle - can be numpy array
185#       which (if either) of these is "right"?
186#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
187#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
188    return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
189        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
190   
191def Oblique(ObCoeff,Tth):
192    if ObCoeff:
193        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
194    else:
195        return 1.0
196               
197def Ruland(RulCoff,wave,Q,Compton):
198    C = 2.9978e8
199    D = 1.5e-3
200    hmc = 0.024262734687
201    sinth2 = (Q*wave/(4.0*np.pi))**2
202    dlam = (wave**2)*Compton*Q/C
203    dlam_c = 2.0*hmc*sinth2-D*wave**2
204    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
205   
206def LorchWeight(Q):
207    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
208           
209def GetAsfMean(ElList,Sthl2):
210#   Calculate various scattering factor terms for PDF calcs
211#   ElList: element dictionary contains scattering factor coefficients, etc.
212#   Sthl2: numpy array of sin theta/lambda squared values
213#   returns: mean(f^2), mean(f)^2, mean(compton)
214    sumNoAtoms = 0.0
215    FF = np.zeros_like(Sthl2)
216    FF2 = np.zeros_like(Sthl2)
217    CF = np.zeros_like(Sthl2)
218    for El in ElList:
219        sumNoAtoms += ElList[El]['FormulaNo']
220    for El in ElList:
221        el = ElList[El]
222        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
223        cf = G2elem.ComptonFac(el,Sthl2)
224        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
225        FF2 += ff2*el['FormulaNo']/sumNoAtoms
226        CF += cf*el['FormulaNo']/sumNoAtoms
227    return FF2,FF**2,CF
228   
229def GetNumDensity(ElList,Vol):
230    sumNoAtoms = 0.0
231    for El in ElList:
232        sumNoAtoms += ElList[El]['FormulaNo']
233    return sumNoAtoms/Vol
234           
235def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
236    NS = 500
237    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
238    Cgama = np.cos(Gama)[:,np.newaxis]
239    Sgama = np.sin(Gama)[:,np.newaxis]
240    CSgama = 1.0/Sgama
241    Delt = Gama[1]-Gama[0]
242    emc = 7.94e-26
243    Navo = 6.023e23
244    Cth = npcosd(Tth/2.0)
245    CTth = npcosd(Tth)
246    Sth = npcosd(Tth/2.0)
247    STth = npsind(Tth)
248    CSth = 1.0/Sth
249    CSTth = 1.0/STth
250    SCth = 1.0/Cth
251    SCTth = 1.0/CTth
252    if 'Bragg' in Geometry:
253        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
254        Y1 = np.pi
255        Y2 = np.pi/2.0
256        Y3 = 3.*np.pi/8. #3pi/4?
257        W = 1.0/(Sth+np.fabs(Sgama))
258        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
259            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
260        Fac0 = Sth**2*Sgama**2
261        X = Fac0+(Fac0+CTth)**2/2
262        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
263        Z = Cgama**4*Cth**4/2.0
264        E = 2.0*(1.0-a)/(b*Cgama/Cth)
265        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
266        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
267        T1 = np.pi/np.sqrt(F1**2-1.0)
268        T2 = np.pi/np.sqrt(F2**2-1.0)
269        Y4 = T1+T2
270        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
271        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
272        Y7 = (T2-T1)/(F1-F2)
273        YT = F2**2*T2-F1**2*T1
274        Y8 = Y1+YT/(F1-F2)
275        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
276        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
277       
278        Q = np.sum(W*M,axis=0)
279        return Q*G*NS/(NS-1.)
280#
281#      cos2th=2.0d*costh^2 - 1.0d
282#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
283#      Y1=3.1415926d
284#      Y2=Y1*0.5d
285#      Y3=Y2*0.75d
286#      for i=1,num_steps-1 do begin
287#         cosgama=double(cos(gama[i]))
288#         singama=double(sin(gama[i]))
289#         cscgama=1.0d / singama
290#
291#         W=1.0d/(sinth+abs(singama))
292#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
293#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
294#
295#         factor0=sinth^2*singama^2
296#         X=factor0+(factor0+cos2th)^2/2.0d
297#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
298#         Z=cosgama^4/2.0d*costh^4
299#         E=2.0d*(1.0-a)/b/cosgama/costh
300#
301#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
302#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
303#
304#         T1=3.14159/sqrt(F1^2-1.0d)
305#         T2=3.14159/sqrt(F2^2-1.0d)
306#         Y4=T1+T2
307#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
308#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
309#         Y7=(T2-T1)/(F1-F2)
310#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
311#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
312#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
313#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
314#
315#         Q=Q+W*M
316#
317#      endfor
318#      Q=double(num_steps)/Double(num_steps-1)*Q*G
319#      end
320    elif 'Cylinder' in Geometry:
321        Q = 0.
322    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
323        EMA = np.exp(-MuT*(1.0-SCTth))
324        Fac1 = (1.-EMA)/(1.0-SCTth)
325        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
326        Fac0 = Cgama/(1-Sgama*SCTth)
327        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
328        Fac0 = Cgama/(1.0+Sgama*SCTth)
329        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
330        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
331        Y = Sgama**3*Cgama*STth*CTth
332        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
333        Fac2 = 1.0/(b*Cgama*STth)
334        U = 2.0*(1.0-a)*Fac2
335        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
336        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)
337        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
338        Y = -Y
339        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)
340        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
341        return Q*G*NS/(NS-1.)
342    elif 'Tilting' in Geometry:
343        EMA = np.exp(-MuT*SCth)
344        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
345#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
346#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
347#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
348#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
349        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
350#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
351        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
352#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
353        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
354        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
355        Z = 0.5*Cgama**4*Sth**4
356#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
357#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
358#
359#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
360#
361#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
362#
363        AlP = 2.0+b*(1.0-Cth*Sgama)
364        AlM = 2.0+b*(1.0+Cth*Sgama)
365#          alphaplus = 2.0 + b*(1.0 - costh*singama)
366#          alphaminus = 2.0 + b*(1.0 + costh*singama)
367        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
368        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
369#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
370#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
371        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)* \
372            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
373            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
374#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
375#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
376#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
377#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
378#                   cosgama*cosgama/b/b*sinth*sinth))
379        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)* \
380            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
381            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
382#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
383#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
384#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
385#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
386#                   cosgama*cosgama/b/b*sinth*sinth))
387        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
388        return Q*G*NS/(NS-1.)
389#       expmutsecth = Exp(-1.0*mut*secth)
390#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
391#       for i=1, num_steps-1 do begin
392#          cosgama=double(cos(gama[i]))
393#          singama=double(sin(gama[i]))
394#          cscgama=1.0d/singama
395#
396#
397#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
398#
399#
400#
401#
402#    ;               print, a,b
403#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
404#          Q=Q+ Wplus*Mplus + Wminus*Mminus
405#      endfor
406#      Q=double(num_steps)/double(num_steps-1)*Q*G
407#   ;   print, min(q), max(q)
408#      end
409
410def MultiScattering(Geometry,ElList,Tth):
411    BN = BD = 0.0
412    Amu = 0.0
413    for El in ElList:
414        el = ElList[El]
415        BN += el['Z']*el['FormulaNo']
416        BD += el['FormulaNo']
417        Amu += el['FormulaNo']*el['mu']
418       
419#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
420
421np.seterr(divide='ignore')
422
423# Voigt function = convolution (norm,cauchy)
424
425class voigt_gen(st.rv_continuous):
426    '''
427    Voigt function
428    Parameters
429    -----------------------------------------
430    '''
431    def _pdf(self,x,s,g):
432        z = np.empty([2,len(x)])
433        z[0] = st.norm.pdf(x,scale=s)
434        z[1] = st.cauchy.pdf(x,scale=g)
435        Z = fft(z)
436        return fftshift(ifft(Z.prod(axis=0))).real
437#    def _cdf(self, x):
438#    def _ppf(self, q):
439#    def _sf(self, x):
440#    def _isf(self, q):
441#    def _stats(self):
442#    def _entropy(self):
443
444voigt = voigt_gen(name='voigt',shapes='ts,g')
445       
446# Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
447
448class fcjde_gen(st.rv_continuous):
449    """
450    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
451    Ref: J. Appl. Cryst. (1994) 27, 892-900.
452    Parameters
453    -----------------------------------------
454    x: array -1 to 1
455    t: 2-theta position of peak
456    s: sum(S/L,H/L); S: sample height, H: detector opening,
457        L: sample to detector opening distance
458    dx: 2-theta step size in deg
459    Result for fcj.pdf
460    -----------------------------------------
461    if x < t & s = S/L+H/L:
462        fcj.pdf = [1/sqrt({cos(x)**2/cos(t)**2}-1) - 1/s]/cos(x)   
463    if x >= t:
464        fcj.pdf = 0   
465    """
466    def _pdf(self,x,t,s,dx):
467#        T = np.where(t<=90.,dx*x+t,180.-dx*x-t)
468        T = dx*x+t
469        ax = npcosd(T)**2
470        bx = npcosd(t)**2
471        bx = np.where(ax>bx,bx,ax)
472        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
473        fx = np.where(fx > 0.,fx,0.0)
474        return fx
475#    def _cdf(self, x):
476#    def _ppf(self, q):
477#    def _sf(self, x):
478#    def _isf(self, q):
479#    def _stats(self):
480#    def _entropy(self):
481       
482fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
483               
484def getBackground(parmDict,bakType,xdata):
485    if bakType == 'chebyschev':
486        yb = np.zeros_like(xdata)
487        iBak = 0
488        while True:
489            key = 'Back'+str(iBak)
490            try:
491                yb += parmDict[key]*(xdata-xdata[0])**iBak
492                iBak += 1
493            except KeyError:
494                return yb
495               
496def getPeakProfile(parmDict,xdata,varyList,bakType):
497   
498    def calcPeakFFT(x,fxns,widths,pos,args):
499        dx = x[1]-x[0]
500        z = np.empty([len(fxns),len(x)])
501        for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
502            z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
503        Z = fft(z)
504        if len(fxns)%2:
505            return ifft(Z.prod(axis=0)).real
506        else:
507            return fftshift(ifft(Z.prod(axis=0))).real
508                       
509    def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
510       
511        DX = xdata[1]-xdata[0]
512        widths,fmin,fmax = getWidths(pos,sig,gam,shl)
513        x = np.linspace(pos-fmin,pos+fmin,1024)
514        if pos > 90:
515            fmin,fmax = [fmax,fmin]         
516        dx = x[1]-x[0]
517        arg = [pos,shl/10.,dx,]
518        Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
519        if len(np.nonzero(Df)[0])>5:
520            fxns = [st.norm,st.cauchy,fcjde]
521            args = [[],[],arg]
522        else:
523            fxns = [st.norm,st.cauchy]
524            args = [[],[]]
525        Df = calcPeakFFT(x,fxns,widths,pos,args)
526        Df /= np.sum(Df)
527        Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
528        return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
529                   
530    yb = getBackground(parmDict,bakType,xdata)
531    yc = np.zeros_like(yb)
532    U = parmDict['U']
533    V = parmDict['V']
534    W = parmDict['W']
535    X = parmDict['X']
536    Y = parmDict['Y']
537    shl = max(parmDict['SH/L'],0.001)
538    Ka2 = False
539    if 'Lam1' in parmDict.keys():
540        Ka2 = True
541        lamRatio = parmDict['Lam2']/parmDict['Lam1']
542        kRatio = parmDict['I(L2)/I(L1)']
543    iPeak = 0
544    while True:
545        try:
546            pos = parmDict['pos'+str(iPeak)]
547            intens = parmDict['int'+str(iPeak)]
548            sigName = 'sig'+str(iPeak)
549            if sigName in varyList:
550                sig = parmDict[sigName]
551            else:
552                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
553            sig = max(sig,0.001)          #avoid neg sigma
554            gamName = 'gam'+str(iPeak)
555            if gamName in varyList:
556                gam = parmDict[gamName]
557            else:
558                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
559            gam = max(gam,0.01)             #avoid neg gamma
560            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
561            if pos > 90:
562                fmin,fmax = [fmax,fmin]         
563            iBeg = np.searchsorted(xdata,pos-fmin)
564            iFin = np.searchsorted(xdata,pos+fmax)
565            if not iBeg+iFin:       #peak below low limit
566                iPeak += 1
567                continue
568            elif not iBeg-iFin:     #peak above high limit
569                return yb+yc
570            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
571            if Ka2:
572                pos2 = 2.0*asind(lamRatio*sind(pos/2.0))
573                Wd,fmin,fmax = getWidths(pos2,sig,gam,shl)
574                if pos > 90:
575                    fmin,fmax = [fmax,fmin]         
576                iBeg = np.searchsorted(xdata,pos2-fmin)
577                iFin = np.searchsorted(xdata,pos2+fmax)
578                yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
579            iPeak += 1
580        except KeyError:        #no more peaks to process
581            return yb+yc
582       
583def getWidths(pos,sig,gam,shl):
584    if shl:
585        widths = [np.sqrt(sig)/100.,gam/200.,.1]
586    else:
587        widths = [np.sqrt(sig)/100.,gam/200.]
588    fwhm = 2.355*widths[0]+2.*widths[1]
589    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
590    fmax = 15.0*fwhm
591    return widths,fmin,fmax
592               
593def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
594   
595    def SetBackgroundParms(Background):
596        bakType,bakFlag = Background[:2]
597        backVals = Background[3:]
598        backNames = ['Back'+str(i) for i in range(len(backVals))]
599        if bakFlag: #returns backNames as varyList = backNames
600            return bakType,dict(zip(backNames,backVals)),backNames
601        else:       #no background varied; varyList = []
602            return bakType,dict(zip(backNames,backVals)),[]
603       
604    def GetBackgroundParms(parmList,Background):
605        iBak = 0
606        while True:
607            try:
608                bakName = 'Back'+str(iBak)
609                Background[iBak+3] = parmList[bakName]
610                iBak += 1
611            except KeyError:
612                break
613               
614    def BackgroundPrint(Background,sigDict):
615        if Background[1]:
616            print 'Background coefficients for',Background[0],'function'
617            ptfmt = "%12.5f"
618            ptstr =  'values:'
619            sigstr = 'esds  :'
620            for i,back in enumerate(Background[3:]):
621                ptstr += ptfmt % (back)
622                sigstr += ptfmt % (sigDict['Back'+str(i)])
623            print ptstr
624            print sigstr
625        else:
626            print 'Background not refined'
627           
628    def SetInstParms(Inst):
629        insVals,insFlags,insNames = Inst[1:4]
630        dataType = insVals[0]
631        insVary = []
632        for i,flag in enumerate(insFlags):
633            if flag:
634                insVary.append(insNames[i])
635        instDict = dict(zip(insNames,insVals))
636        instDict['X'] = max(instDict['X'],0.1)
637        instDict['Y'] = max(instDict['Y'],0.1)
638        instDict['SH/L'] = max(instDict['SH/L'],0.001)
639        return dataType,instDict,insVary
640       
641    def GetInstParms(parmDict,Inst,varyList,Peaks):
642        instNames = Inst[3]
643        for i,name in enumerate(instNames):
644            Inst[1][i] = parmDict[name]
645        iPeak = 0
646        while True:
647            try:
648                sigName = 'sig'+str(iPeak)
649                pos = parmDict['pos'+str(iPeak)]
650                if sigName not in varyList:
651                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
652                gamName = 'gam'+str(iPeak)
653                if gamName not in varyList:
654                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
655                iPeak += 1
656            except KeyError:
657                break
658       
659    def InstPrint(Inst,sigDict):
660        print 'Instrument Parameters:'
661        ptfmt = "%12.6f"
662        ptlbls = 'names :'
663        ptstr =  'values:'
664        sigstr = 'esds  :'
665        instNames = Inst[3][1:]
666        for i,name in enumerate(instNames):
667            ptlbls += "%s" % (name.center(12))
668            ptstr += ptfmt % (Inst[1][i+1])
669            if name in sigDict:
670                sigstr += ptfmt % (sigDict[name])
671            else:
672                sigstr += 12*' '
673        print ptlbls
674        print ptstr
675        print sigstr
676
677    def setPeaksParms(Peaks):
678        peakNames = []
679        peakVary = []
680        peakVals = []
681        names = ['pos','int','sig','gam']
682        for i,peak in enumerate(Peaks):
683            for j in range(4):
684                peakVals.append(peak[2*j])
685                parName = names[j]+str(i)
686                peakNames.append(parName)
687                if peak[2*j+1]:
688                    peakVary.append(parName)
689        return dict(zip(peakNames,peakVals)),peakVary
690               
691    def GetPeaksParms(parmDict,Peaks,varyList):
692        names = ['pos','int','sig','gam']
693        for i,peak in enumerate(Peaks):
694            for j in range(4):
695                pos = parmDict['pos'+str(i)]
696                parName = names[j]+str(i)
697                if parName in varyList:
698                    peak[2*j] = parmDict[parName]
699                elif 'sig' in parName:
700                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
701                elif 'gam' in parName:
702                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
703                       
704    def PeaksPrint(parmDict,sigDict,varyList):
705        print 'Peak coefficients:'
706        names = ['pos','int','sig','gam']
707        head = 15*' '
708        for name in names:
709            head += name.center(12)+'esd'.center(12)
710        print head
711        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
712        for i,peak in enumerate(Peaks):
713            ptstr =  ':'
714            for j in range(4):
715                name = names[j]
716                parName = name+str(i)
717                ptstr += ptfmt[name] % (parmDict[parName])
718                if parName in varyList:
719#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
720                    ptstr += ptfmt[name] % (sigDict[parName])
721                else:
722#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
723                    ptstr += 12*' '
724            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
725           
726    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
727        parmdict.update(zip(varylist,values))
728        M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
729        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
730        if dlg:
731            dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))
732        return M
733       
734    x,y,w,yc,yb,yd = data               #these are numpy arrays!
735    xBeg = np.searchsorted(x,Limits[0])
736    xFin = np.searchsorted(x,Limits[1])
737    bakType,bakDict,bakVary = SetBackgroundParms(Background)
738    dataType,insDict,insVary = SetInstParms(Inst)
739    peakDict,peakVary = setPeaksParms(Peaks)
740    parmDict = {}
741    parmDict.update(bakDict)
742    parmDict.update(insDict)
743    parmDict.update(peakDict)
744    varyList = bakVary+insVary+peakVary
745    while True:
746        begin = time.time()
747        values =  np.array(ValuesOut(parmDict, varyList))
748        if FitPgm == 'LSQ':
749            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
750            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME)
751            screenSize = wx.ClientDisplayRect()
752            Size = dlg.GetSize()
753            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
754            try:
755                result = so.leastsq(errPeakProfile,values,
756                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
757            finally:
758                dlg.Destroy()
759            runtime = time.time()-begin   
760            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
761            print "%s%8.3f%s " % ('fitpeak time =',runtime,'s')
762            ValuesIn(parmDict, varyList, result[0])
763            chisq = np.sum(errPeakProfile(result[0],x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,0)**2)
764            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
765            GOF = chisq/(xFin-xBeg-len(varyList))
766            print "%s%7.2f%s%12.6g%s%6.2f" % ('Rwp = ',Rwp,'%, chi**2 = ',chisq,' reduced chi**2 = ',GOF)
767            try:
768                sig = np.sqrt(np.diag(result[1])*chisq)
769                break                   #refinement succeeded - finish up!
770            except ValueError:          #result[1] is None on singular matrix
771                print 'Refinement failed - singular matrix'
772                Ipvt = result[2]['ipvt']
773                for i,ipvt in enumerate(Ipvt):
774                    if not np.sum(result[2]['fjac'],axis=1)[i]:
775                        print 'Removing parameter: ',varyList[ipvt-1]
776                        del(varyList[ipvt-1])
777                        break
778        elif FitPgm == 'BFGS':
779            print 'Other program here'
780            return
781       
782    sigDict = dict(zip(varyList,sig))
783    yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
784    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
785    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
786    GetBackgroundParms(parmDict,Background)
787    BackgroundPrint(Background,sigDict)
788    GetInstParms(parmDict,Inst,varyList,Peaks)
789    InstPrint(Inst,sigDict)
790    GetPeaksParms(parmDict,Peaks,varyList)   
791    PeaksPrint(parmDict,sigDict,varyList)
792   
793def CalcPDF(data,inst,xydata):
794    auxPlot = []
795    import copy
796    import scipy.fftpack as ft
797    #subtract backgrounds - if any
798    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
799    if data['Sample Bkg.']['Name']:
800        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
801            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
802    if data['Container']['Name']:
803        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
804        if data['Container Bkg.']['Name']:
805            xycontainer += (xydata['Container Bkg.'][1][1]+
806                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
807        xydata['IofQ'][1][1] += xycontainer
808    #get element data & absorption coeff.
809    ElList = data['ElList']
810    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
811    #Apply angle dependent corrections
812    Tth = xydata['Sample'][1][0]
813    dt = (Tth[1]-Tth[0])
814    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
815    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
816    if data['DetType'] == 'Image plate':
817        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
818    XY = xydata['IofQ'][1]   
819    #convert to Q
820    hc = 12.397639
821    if 'Lam' in inst:
822        wave = inst['Lam']
823    else:
824        wave = inst['Lam1']
825    keV = hc/wave
826    minQ = npT2q(Tth[0],wave)
827    maxQ = npT2q(Tth[-1],wave)   
828    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
829    dq = Qpoints[1]-Qpoints[0]
830    XY[0] = npT2q(XY[0],wave)   
831#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
832    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
833    Qdata = T(Qpoints)
834   
835    qLimits = data['QScaleLim']
836    minQ = np.searchsorted(Qpoints,qLimits[0])
837    maxQ = np.searchsorted(Qpoints,qLimits[1])
838    newdata = []
839    xydata['IofQ'][1][0] = Qpoints
840    xydata['IofQ'][1][1] = Qdata
841    for item in xydata['IofQ'][1]:
842        newdata.append(item[:maxQ])
843    xydata['IofQ'][1] = newdata
844   
845
846    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
847    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
848    Q = xydata['SofQ'][1][0]
849    ruland = Ruland(data['Ruland'],wave,Q,CF)
850#    auxPlot.append([Q,ruland,'Ruland'])     
851    CF *= ruland
852#    auxPlot.append([Q,CF,'CF-Corr'])
853    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
854    xydata['SofQ'][1][1] *= scale
855    xydata['SofQ'][1][1] -= CF
856    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
857    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
858    xydata['SofQ'][1][1] *= scale
859   
860    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
861    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
862    if data['Lorch']:
863        xydata['FofQ'][1][1] *= LorchWeight(Q)
864   
865    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
866    nR = len(xydata['GofR'][1][1])
867    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
868    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
869   
870       
871    return auxPlot
872       
873#testing data
874import plot
875NeedTestData = True
876def TestData():
877#    global NeedTestData
878    NeedTestData = False
879    global bakType
880    bakType = 'chebyschev'
881    global xdata
882    xdata = np.linspace(4.0,40.0,36000)
883    global parmDict0
884    parmDict0 = {
885        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
886        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
887        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
888        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
889        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
890        'Back0':5.384,'Back1':-0.015,'Back2':.004,
891        }
892    global parmDict1
893    parmDict1 = {
894        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
895        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
896        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
897        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
898        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
899        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
900        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
901        'Back0':36.897,'Back1':-0.508,'Back2':.006,
902        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
903        }
904    global varyList
905    varyList = []
906
907def test0():
908    if NeedTestData: TestData()
909    msg = 'test '
910    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
911    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
912    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
913    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
914    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
915    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
916    time0 = time.time()
917    for i in range(100):
918        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
919    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
920   
921
922if __name__ == '__main__':
923    global plotter
924    plotter = plot.PlotNotebook()
925    test0()
926    print "OK"
927    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.