source: trunk/GSASIIpwd.py @ 301

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

Add cancel to progress bar & associated things to make it work
File size: 36.5 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    T = x*dx+t
462    s = S/L+H/L
463    if x < 0:
464        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
465    if x >= 0:
466        fcj.pdf = 0   
467    """
468    def _pdf(self,x,t,s,dx):
469        T = dx*x+t
470        ax = npcosd(T)**2
471        bx = npcosd(t)**2
472        bx = np.where(ax>bx,bx,ax)
473        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
474        fx = np.where(fx > 0.,fx,0.0)
475        return fx
476#    def _cdf(self, x):
477#    def _ppf(self, q):
478#    def _sf(self, x):
479#    def _isf(self, q):
480#    def _stats(self):
481#    def _entropy(self):
482       
483fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
484               
485def getBackground(parmDict,bakType,xdata):
486    if bakType == 'chebyschev':
487        yb = np.zeros_like(xdata)
488        iBak = 0
489        while True:
490            key = 'Back'+str(iBak)
491            try:
492                yb += parmDict[key]*(xdata-xdata[0])**iBak
493                iBak += 1
494            except KeyError:
495                return yb
496               
497def getPeakProfile(parmDict,xdata,varyList,bakType):
498   
499    def calcPeakFFT(x,fxns,widths,pos,args):
500        dx = x[1]-x[0]
501        z = np.empty([len(fxns),len(x)])
502        for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
503            z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
504        Z = fft(z)
505        if len(fxns)%2:
506            return ifft(Z.prod(axis=0)).real
507        else:
508            return fftshift(ifft(Z.prod(axis=0))).real
509                       
510    def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
511       
512        DX = xdata[1]-xdata[0]
513        widths,fmin,fmax = getWidths(pos,sig,gam,shl)
514        x = np.linspace(pos-fmin,pos+fmin,1024)
515        if pos > 90:
516            fmin,fmax = [fmax,fmin]         
517        dx = x[1]-x[0]
518        arg = [pos,shl/10.,dx,]
519        Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
520        if len(np.nonzero(Df)[0])>5:
521            fxns = [st.norm,st.cauchy,fcjde]
522            args = [[],[],arg]
523        else:
524            fxns = [st.norm,st.cauchy]
525            args = [[],[]]
526        Df = calcPeakFFT(x,fxns,widths,pos,args)
527        Df /= np.sum(Df)
528        Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
529        return intens*Df(xdata)*DX/dx       #*10 to get close to old fxn???
530                   
531    yb = getBackground(parmDict,bakType,xdata)
532    yc = np.zeros_like(yb)
533    U = parmDict['U']
534    V = parmDict['V']
535    W = parmDict['W']
536    X = parmDict['X']
537    Y = parmDict['Y']
538    shl = max(parmDict['SH/L'],0.001)
539    Ka2 = False
540    if 'Lam1' in parmDict.keys():
541        Ka2 = True
542        lamRatio = parmDict['Lam2']/parmDict['Lam1']
543        kRatio = parmDict['I(L2)/I(L1)']
544    iPeak = 0
545    while True:
546        try:
547            pos = parmDict['pos'+str(iPeak)]
548            intens = parmDict['int'+str(iPeak)]
549            sigName = 'sig'+str(iPeak)
550            if sigName in varyList:
551                sig = parmDict[sigName]
552            else:
553                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
554            sig = max(sig,0.001)          #avoid neg sigma
555            gamName = 'gam'+str(iPeak)
556            if gamName in varyList:
557                gam = parmDict[gamName]
558            else:
559                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
560            gam = max(gam,0.01)             #avoid neg gamma
561            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
562            if pos > 90:
563                fmin,fmax = [fmax,fmin]         
564            iBeg = np.searchsorted(xdata,pos-fmin)
565            iFin = np.searchsorted(xdata,pos+fmax)
566            if not iBeg+iFin:       #peak below low limit
567                iPeak += 1
568                continue
569            elif not iBeg-iFin:     #peak above high limit
570                return yb+yc
571            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
572            if Ka2:
573                pos2 = 2.0*asind(lamRatio*sind(pos/2.0))
574                Wd,fmin,fmax = getWidths(pos2,sig,gam,shl)
575                if pos > 90:
576                    fmin,fmax = [fmax,fmin]         
577                iBeg = np.searchsorted(xdata,pos2-fmin)
578                iFin = np.searchsorted(xdata,pos2+fmax)
579                yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
580            iPeak += 1
581        except KeyError:        #no more peaks to process
582            return yb+yc
583       
584def getWidths(pos,sig,gam,shl):
585    if shl:
586        widths = [np.sqrt(sig)/100.,gam/200.,.1]
587    else:
588        widths = [np.sqrt(sig)/100.,gam/200.]
589    fwhm = 2.355*widths[0]+2.*widths[1]
590    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
591    fmax = 15.0*fwhm
592    return widths,fmin,fmax
593               
594def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data):
595   
596    def SetBackgroundParms(Background):
597        bakType,bakFlag = Background[:2]
598        backVals = Background[3:]
599        backNames = ['Back'+str(i) for i in range(len(backVals))]
600        if bakFlag: #returns backNames as varyList = backNames
601            return bakType,dict(zip(backNames,backVals)),backNames
602        else:       #no background varied; varyList = []
603            return bakType,dict(zip(backNames,backVals)),[]
604       
605    def GetBackgroundParms(parmList,Background):
606        iBak = 0
607        while True:
608            try:
609                bakName = 'Back'+str(iBak)
610                Background[iBak+3] = parmList[bakName]
611                iBak += 1
612            except KeyError:
613                break
614               
615    def BackgroundPrint(Background,sigDict):
616        if Background[1]:
617            print 'Background coefficients for',Background[0],'function'
618            ptfmt = "%12.5f"
619            ptstr =  'values:'
620            sigstr = 'esds  :'
621            for i,back in enumerate(Background[3:]):
622                ptstr += ptfmt % (back)
623                sigstr += ptfmt % (sigDict['Back'+str(i)])
624            print ptstr
625            print sigstr
626        else:
627            print 'Background not refined'
628           
629    def SetInstParms(Inst):
630        insVals,insFlags,insNames = Inst[1:4]
631        dataType = insVals[0]
632        insVary = []
633        for i,flag in enumerate(insFlags):
634            if flag:
635                insVary.append(insNames[i])
636        instDict = dict(zip(insNames,insVals))
637        instDict['X'] = max(instDict['X'],0.1)
638        instDict['Y'] = max(instDict['Y'],0.1)
639        instDict['SH/L'] = max(instDict['SH/L'],0.001)
640        return dataType,instDict,insVary
641       
642    def GetInstParms(parmDict,Inst,varyList,Peaks):
643        instNames = Inst[3]
644        for i,name in enumerate(instNames):
645            Inst[1][i] = parmDict[name]
646        iPeak = 0
647        while True:
648            try:
649                sigName = 'sig'+str(iPeak)
650                pos = parmDict['pos'+str(iPeak)]
651                if sigName not in varyList:
652                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
653                gamName = 'gam'+str(iPeak)
654                if gamName not in varyList:
655                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
656                iPeak += 1
657            except KeyError:
658                break
659       
660    def InstPrint(Inst,sigDict):
661        print 'Instrument Parameters:'
662        ptfmt = "%12.6f"
663        ptlbls = 'names :'
664        ptstr =  'values:'
665        sigstr = 'esds  :'
666        instNames = Inst[3][1:]
667        for i,name in enumerate(instNames):
668            ptlbls += "%s" % (name.center(12))
669            ptstr += ptfmt % (Inst[1][i+1])
670            if name in sigDict:
671                sigstr += ptfmt % (sigDict[name])
672            else:
673                sigstr += 12*' '
674        print ptlbls
675        print ptstr
676        print sigstr
677
678    def setPeaksParms(Peaks):
679        peakNames = []
680        peakVary = []
681        peakVals = []
682        names = ['pos','int','sig','gam']
683        for i,peak in enumerate(Peaks):
684            for j in range(4):
685                peakVals.append(peak[2*j])
686                parName = names[j]+str(i)
687                peakNames.append(parName)
688                if peak[2*j+1]:
689                    peakVary.append(parName)
690        return dict(zip(peakNames,peakVals)),peakVary
691               
692    def GetPeaksParms(parmDict,Peaks,varyList):
693        names = ['pos','int','sig','gam']
694        for i,peak in enumerate(Peaks):
695            for j in range(4):
696                pos = parmDict['pos'+str(i)]
697                parName = names[j]+str(i)
698                if parName in varyList:
699                    peak[2*j] = parmDict[parName]
700                elif 'sig' in parName:
701                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
702                elif 'gam' in parName:
703                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
704                       
705    def PeaksPrint(parmDict,sigDict,varyList):
706        print 'Peak coefficients:'
707        names = ['pos','int','sig','gam']
708        head = 15*' '
709        for name in names:
710            head += name.center(12)+'esd'.center(12)
711        print head
712        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
713        for i,peak in enumerate(Peaks):
714            ptstr =  ':'
715            for j in range(4):
716                name = names[j]
717                parName = name+str(i)
718                ptstr += ptfmt[name] % (parmDict[parName])
719                if parName in varyList:
720#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
721                    ptstr += ptfmt[name] % (sigDict[parName])
722                else:
723#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
724                    ptstr += 12*' '
725            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
726           
727    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
728        parmdict.update(zip(varylist,values))
729        M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
730        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
731        if dlg:
732            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
733            if not GoOn:
734                return -M           #abort!!
735        return M
736       
737    x,y,w,yc,yb,yd = data               #these are numpy arrays!
738    xBeg = np.searchsorted(x,Limits[0])
739    xFin = np.searchsorted(x,Limits[1])
740    bakType,bakDict,bakVary = SetBackgroundParms(Background)
741    dataType,insDict,insVary = SetInstParms(Inst)
742    peakDict,peakVary = setPeaksParms(Peaks)
743    parmDict = {}
744    parmDict.update(bakDict)
745    parmDict.update(insDict)
746    parmDict.update(peakDict)
747    varyList = bakVary+insVary+peakVary
748    while True:
749        begin = time.time()
750        values =  np.array(ValuesOut(parmDict, varyList))
751        if FitPgm == 'LSQ':
752            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
753            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
754            screenSize = wx.ClientDisplayRect()
755            Size = dlg.GetSize()
756            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
757            try:
758                result = so.leastsq(errPeakProfile,values,
759                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
760            finally:
761                dlg.Destroy()
762            runtime = time.time()-begin   
763            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
764            print "%s%8.3f%s " % ('fitpeak time =',runtime,'s')
765            ValuesIn(parmDict, varyList, result[0])
766            chisq = np.sum(errPeakProfile(result[0],x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,0)**2)
767            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
768            GOF = chisq/(xFin-xBeg-len(varyList))
769            print "%s%7.2f%s%12.6g%s%6.2f" % ('Rwp = ',Rwp,'%, chi**2 = ',chisq,' reduced chi**2 = ',GOF)
770            try:
771                sig = np.sqrt(np.diag(result[1])*chisq)
772                if np.any(np.isnan(sig)):
773                    print '*** Least squares aborted - some invalid esds possible ***'
774                break                   #refinement succeeded - finish up!
775            except ValueError:          #result[1] is None on singular matrix
776                print 'Refinement failed - singular matrix'
777                Ipvt = result[2]['ipvt']
778                for i,ipvt in enumerate(Ipvt):
779                    if not np.sum(result[2]['fjac'],axis=1)[i]:
780                        print 'Removing parameter: ',varyList[ipvt-1]
781                        del(varyList[ipvt-1])
782                        break
783        elif FitPgm == 'BFGS':
784            print 'Other program here'
785            return
786       
787    sigDict = dict(zip(varyList,sig))
788    yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
789    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
790    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
791    GetBackgroundParms(parmDict,Background)
792    BackgroundPrint(Background,sigDict)
793    GetInstParms(parmDict,Inst,varyList,Peaks)
794    InstPrint(Inst,sigDict)
795    GetPeaksParms(parmDict,Peaks,varyList)   
796    PeaksPrint(parmDict,sigDict,varyList)
797   
798def CalcPDF(data,inst,xydata):
799    auxPlot = []
800    import copy
801    import scipy.fftpack as ft
802    #subtract backgrounds - if any
803    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
804    if data['Sample Bkg.']['Name']:
805        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
806            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
807    if data['Container']['Name']:
808        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
809        if data['Container Bkg.']['Name']:
810            xycontainer += (xydata['Container Bkg.'][1][1]+
811                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
812        xydata['IofQ'][1][1] += xycontainer
813    #get element data & absorption coeff.
814    ElList = data['ElList']
815    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
816    #Apply angle dependent corrections
817    Tth = xydata['Sample'][1][0]
818    dt = (Tth[1]-Tth[0])
819    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
820    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
821    if data['DetType'] == 'Image plate':
822        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
823    XY = xydata['IofQ'][1]   
824    #convert to Q
825    hc = 12.397639
826    if 'Lam' in inst:
827        wave = inst['Lam']
828    else:
829        wave = inst['Lam1']
830    keV = hc/wave
831    minQ = npT2q(Tth[0],wave)
832    maxQ = npT2q(Tth[-1],wave)   
833    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
834    dq = Qpoints[1]-Qpoints[0]
835    XY[0] = npT2q(XY[0],wave)   
836#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
837    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
838    Qdata = T(Qpoints)
839   
840    qLimits = data['QScaleLim']
841    minQ = np.searchsorted(Qpoints,qLimits[0])
842    maxQ = np.searchsorted(Qpoints,qLimits[1])
843    newdata = []
844    xydata['IofQ'][1][0] = Qpoints
845    xydata['IofQ'][1][1] = Qdata
846    for item in xydata['IofQ'][1]:
847        newdata.append(item[:maxQ])
848    xydata['IofQ'][1] = newdata
849   
850
851    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
852    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
853    Q = xydata['SofQ'][1][0]
854    ruland = Ruland(data['Ruland'],wave,Q,CF)
855#    auxPlot.append([Q,ruland,'Ruland'])     
856    CF *= ruland
857#    auxPlot.append([Q,CF,'CF-Corr'])
858    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
859    xydata['SofQ'][1][1] *= scale
860    xydata['SofQ'][1][1] -= CF
861    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
862    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
863    xydata['SofQ'][1][1] *= scale
864   
865    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
866    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
867    if data['Lorch']:
868        xydata['FofQ'][1][1] *= LorchWeight(Q)
869   
870    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
871    nR = len(xydata['GofR'][1][1])
872    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
873    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
874   
875       
876    return auxPlot
877       
878#testing data
879import plot
880NeedTestData = True
881def TestData():
882#    global NeedTestData
883    NeedTestData = False
884    global bakType
885    bakType = 'chebyschev'
886    global xdata
887    xdata = np.linspace(4.0,40.0,36000)
888    global parmDict0
889    parmDict0 = {
890        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
891        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
892        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
893        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
894        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
895        'Back0':5.384,'Back1':-0.015,'Back2':.004,
896        }
897    global parmDict1
898    parmDict1 = {
899        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
900        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
901        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
902        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
903        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
904        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
905        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
906        'Back0':36.897,'Back1':-0.508,'Back2':.006,
907        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
908        }
909    global varyList
910    varyList = []
911
912def test0():
913    if NeedTestData: TestData()
914    msg = 'test '
915    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
916    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
917    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
918    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
919    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
920    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
921    time0 = time.time()
922    for i in range(100):
923        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
924    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
925   
926
927if __name__ == '__main__':
928    global plotter
929    plotter = plot.PlotNotebook()
930    test0()
931    print "OK"
932    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.