source: trunk/GSASIIpwd.py @ 315

Last change on this file since 315 was 315, checked in by toby, 12 years ago

new scons file for compiling (needs linux & win work); update macbin's; changes to texture to work on Mac; minor cleanup to GSASIIpwd
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
25import GSASIIplot as G2plt
26import GSASIIlattice as G2lat
27import GSASIIElem as G2elem
28import GSASIIgrid as G2gd
29import GSASIIIO as G2IO
30
31# trig functions in degrees
32sind = lambda x: math.sin(x*math.pi/180.)
33asind = lambda x: 180.*math.asin(x)/math.pi
34tand = lambda x: math.tan(x*math.pi/180.)
35atand = lambda x: 180.*math.atan(x)/math.pi
36atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
37cosd = lambda x: math.cos(x*math.pi/180.)
38acosd = lambda x: 180.*math.acos(x)/math.pi
39rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
40#numpy versions
41npsind = lambda x: np.sin(x*np.pi/180.)
42npasind = lambda x: 180.*np.arcsin(x)/math.pi
43npcosd = lambda x: np.cos(x*math.pi/180.)
44npacosd = lambda x: 180.*np.arccos(x)/math.pi
45nptand = lambda x: np.tan(x*math.pi/180.)
46npatand = lambda x: 180.*np.arctan(x)/np.pi
47npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
48npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
49npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
50   
51       
52def factorize(num):
53    ''' Provide prime number factors for integer num
54    Returns dictionary of prime factors (keys) & power for each (data)
55    '''
56    factors = {}
57    orig = num
58
59    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
60    i, sqi = 2, 4
61    while sqi <= num:
62        while not num%i:
63            num /= i
64            factors[i] = factors.get(i, 0) + 1
65
66        sqi += 2*i + 1
67        i += 1
68
69    if num != 1 and num != orig:
70        factors[num] = factors.get(num, 0) + 1
71
72    if factors:
73        return factors
74    else:
75        return {num:1}          #a prime number!
76           
77def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
78    ''' Provide list of optimal data sizes for FFT calculations
79    Input:
80        nmin: minimum data size >= 1
81        nmax: maximum data size > nmin
82        thresh: maximum prime factor allowed
83    Returns:
84        list of data sizes where the maximum prime factor is < thresh
85    ''' 
86    plist = []
87    nmin = max(1,nmin)
88    nmax = max(nmin+1,nmax)
89    for p in range(nmin,nmax):
90        if max(factorize(p).keys()) < thresh:
91            plist.append(p)
92    return plist
93
94def ValuesOut(parmdict, varylist):
95    '''Use before call to leastsq to setup list of values for the parameters
96    in parmdict, as selected by key in varylist'''
97    return [parmdict[key] for key in varylist] 
98   
99def ValuesIn(parmdict, varylist, values):
100    ''' Use after call to leastsq to update the parameter dictionary with
101    values corresponding to keys in varylist'''
102    parmdict.update(zip(varylist,values))
103   
104def Transmission(Geometry,Abs,Diam):
105#Calculate sample transmission
106#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
107#   Abs: absorption coeff in cm-1
108#   Diam: sample thickness/diameter in mm
109    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
110        MuR = Abs*Diam/20.0
111        if MuR <= 3.0:
112            T0 = 16/(3.*math.pi)
113            T1 = -0.045780
114            T2 = -0.02489
115            T3 = 0.003045
116            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
117            if T < -20.:
118                return 2.06e-9
119            else:
120                return math.exp(T)
121        else:
122            T1 = 1.433902
123            T2 = 0.013869+0.337894
124            T3 = 1.933433+1.163198
125            T4 = 0.044365-0.04259
126            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
127            return T/100.
128    elif 'plate' in Geometry:
129        MuR = Abs*Diam/10.
130        return math.exp(-MuR)
131    elif 'Bragg' in Geometry:
132        return 0.0
133
134def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
135#Calculate sample absorption
136#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
137#   Abs: absorption coeff in cm-1
138#   Diam: sample thickness/diameter in mm
139#   Tth: 2-theta scattering angle - can be numpy array
140#   Phi: flat plate tilt angle - future
141#   Psi: flat plate tilt axis - future
142    Sth2 = npsind(Tth/2.0)**2
143    Cth2 = 1.-Sth2
144    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
145        MuR = Abs*Diam/20.0
146        if MuR < 3.0:
147            T0 = 16.0/(3*np.pi)
148            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
149                0.109561*np.sqrt(Sth2)-26.04556
150            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
151                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
152            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
153            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
154            return np.exp(Trns)
155        else:
156            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
157                10.02088*Sth2**3-3.36778*Sth2**4
158            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
159                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
160            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
161                0.13576*np.sqrt(Sth2)+1.163198
162            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
163            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
164            return Trns/100.
165    elif 'Bragg' in Geometry:
166        return 1.0
167    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
168        # and only defined for 2theta < 90
169        MuR = Abs*Diam/10.0
170        T1 = np.exp(-MuR)
171        T2 = np.exp(-MuR/npcosd(Tth))
172        Tb = MuR-MuR/npcosd(Tth)
173        return (T2-T1)/Tb
174    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
175        MuR = Abs*Diam/10.0
176        cth = npcosd(Tth/2.0)
177        return np.exp(-MuR/cth)/cth
178       
179def Polarization(Pola,Tth,Azm=0.0):
180#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
181#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
182#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
183#   Tth: 2-theta scattering angle - can be numpy array
184#       which (if either) of these is "right"?
185#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
186#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
187    return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
188        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
189   
190def Oblique(ObCoeff,Tth):
191    if ObCoeff:
192        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
193    else:
194        return 1.0
195               
196def Ruland(RulCoff,wave,Q,Compton):
197    C = 2.9978e8
198    D = 1.5e-3
199    hmc = 0.024262734687
200    sinth2 = (Q*wave/(4.0*np.pi))**2
201    dlam = (wave**2)*Compton*Q/C
202    dlam_c = 2.0*hmc*sinth2-D*wave**2
203    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
204   
205def LorchWeight(Q):
206    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
207           
208def GetAsfMean(ElList,Sthl2):
209#   Calculate various scattering factor terms for PDF calcs
210#   ElList: element dictionary contains scattering factor coefficients, etc.
211#   Sthl2: numpy array of sin theta/lambda squared values
212#   returns: mean(f^2), mean(f)^2, mean(compton)
213    sumNoAtoms = 0.0
214    FF = np.zeros_like(Sthl2)
215    FF2 = np.zeros_like(Sthl2)
216    CF = np.zeros_like(Sthl2)
217    for El in ElList:
218        sumNoAtoms += ElList[El]['FormulaNo']
219    for El in ElList:
220        el = ElList[El]
221        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
222        cf = G2elem.ComptonFac(el,Sthl2)
223        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
224        FF2 += ff2*el['FormulaNo']/sumNoAtoms
225        CF += cf*el['FormulaNo']/sumNoAtoms
226    return FF2,FF**2,CF
227   
228def GetNumDensity(ElList,Vol):
229    sumNoAtoms = 0.0
230    for El in ElList:
231        sumNoAtoms += ElList[El]['FormulaNo']
232    return sumNoAtoms/Vol
233           
234def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
235    NS = 500
236    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
237    Cgama = np.cos(Gama)[:,np.newaxis]
238    Sgama = np.sin(Gama)[:,np.newaxis]
239    CSgama = 1.0/Sgama
240    Delt = Gama[1]-Gama[0]
241    emc = 7.94e-26
242    Navo = 6.023e23
243    Cth = npcosd(Tth/2.0)
244    CTth = npcosd(Tth)
245    Sth = npcosd(Tth/2.0)
246    STth = npsind(Tth)
247    CSth = 1.0/Sth
248    CSTth = 1.0/STth
249    SCth = 1.0/Cth
250    SCTth = 1.0/CTth
251    if 'Bragg' in Geometry:
252        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
253        Y1 = np.pi
254        Y2 = np.pi/2.0
255        Y3 = 3.*np.pi/8. #3pi/4?
256        W = 1.0/(Sth+np.fabs(Sgama))
257        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
258            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
259        Fac0 = Sth**2*Sgama**2
260        X = Fac0+(Fac0+CTth)**2/2
261        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
262        Z = Cgama**4*Cth**4/2.0
263        E = 2.0*(1.0-a)/(b*Cgama/Cth)
264        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
265        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
266        T1 = np.pi/np.sqrt(F1**2-1.0)
267        T2 = np.pi/np.sqrt(F2**2-1.0)
268        Y4 = T1+T2
269        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
270        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
271        Y7 = (T2-T1)/(F1-F2)
272        YT = F2**2*T2-F1**2*T1
273        Y8 = Y1+YT/(F1-F2)
274        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
275        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
276       
277        Q = np.sum(W*M,axis=0)
278        return Q*G*NS/(NS-1.)
279#
280#      cos2th=2.0d*costh^2 - 1.0d
281#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
282#      Y1=3.1415926d
283#      Y2=Y1*0.5d
284#      Y3=Y2*0.75d
285#      for i=1,num_steps-1 do begin
286#         cosgama=double(cos(gama[i]))
287#         singama=double(sin(gama[i]))
288#         cscgama=1.0d / singama
289#
290#         W=1.0d/(sinth+abs(singama))
291#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
292#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
293#
294#         factor0=sinth^2*singama^2
295#         X=factor0+(factor0+cos2th)^2/2.0d
296#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
297#         Z=cosgama^4/2.0d*costh^4
298#         E=2.0d*(1.0-a)/b/cosgama/costh
299#
300#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
301#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
302#
303#         T1=3.14159/sqrt(F1^2-1.0d)
304#         T2=3.14159/sqrt(F2^2-1.0d)
305#         Y4=T1+T2
306#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
307#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
308#         Y7=(T2-T1)/(F1-F2)
309#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
310#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
311#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
312#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
313#
314#         Q=Q+W*M
315#
316#      endfor
317#      Q=double(num_steps)/Double(num_steps-1)*Q*G
318#      end
319    elif 'Cylinder' in Geometry:
320        Q = 0.
321    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
322        EMA = np.exp(-MuT*(1.0-SCTth))
323        Fac1 = (1.-EMA)/(1.0-SCTth)
324        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
325        Fac0 = Cgama/(1-Sgama*SCTth)
326        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
327        Fac0 = Cgama/(1.0+Sgama*SCTth)
328        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
329        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
330        Y = Sgama**3*Cgama*STth*CTth
331        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
332        Fac2 = 1.0/(b*Cgama*STth)
333        U = 2.0*(1.0-a)*Fac2
334        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
335        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)
336        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
337        Y = -Y
338        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)
339        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
340        return Q*G*NS/(NS-1.)
341    elif 'Tilting' in Geometry:
342        EMA = np.exp(-MuT*SCth)
343        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
344#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
345#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
346#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
347#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
348        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
349#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
350        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
351#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
352        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
353        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
354        Z = 0.5*Cgama**4*Sth**4
355#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
356#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
357#
358#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
359#
360#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
361#
362        AlP = 2.0+b*(1.0-Cth*Sgama)
363        AlM = 2.0+b*(1.0+Cth*Sgama)
364#          alphaplus = 2.0 + b*(1.0 - costh*singama)
365#          alphaminus = 2.0 + b*(1.0 + costh*singama)
366        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
367        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
368#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
369#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
370        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)* \
371            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
372            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
373#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
374#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
375#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
376#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
377#                   cosgama*cosgama/b/b*sinth*sinth))
378        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)* \
379            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
380            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
381#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
382#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
383#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
384#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
385#                   cosgama*cosgama/b/b*sinth*sinth))
386        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
387        return Q*G*NS/(NS-1.)
388#       expmutsecth = Exp(-1.0*mut*secth)
389#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
390#       for i=1, num_steps-1 do begin
391#          cosgama=double(cos(gama[i]))
392#          singama=double(sin(gama[i]))
393#          cscgama=1.0d/singama
394#
395#
396#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
397#
398#
399#
400#
401#    ;               print, a,b
402#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
403#          Q=Q+ Wplus*Mplus + Wminus*Mminus
404#      endfor
405#      Q=double(num_steps)/double(num_steps-1)*Q*G
406#   ;   print, min(q), max(q)
407#      end
408
409def MultiScattering(Geometry,ElList,Tth):
410    BN = BD = 0.0
411    Amu = 0.0
412    for El in ElList:
413        el = ElList[El]
414        BN += el['Z']*el['FormulaNo']
415        BD += el['FormulaNo']
416        Amu += el['FormulaNo']*el['mu']
417       
418#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
419
420np.seterr(divide='ignore')
421
422# Voigt function = convolution (norm,cauchy)
423
424class voigt_gen(st.rv_continuous):
425    '''
426    Voigt function
427    Parameters
428    -----------------------------------------
429    '''
430    def _pdf(self,x,s,g):
431        z = np.empty([2,len(x)])
432        z[0] = st.norm.pdf(x,scale=s)
433        z[1] = st.cauchy.pdf(x,scale=g)
434        Z = fft(z)
435        return fftshift(ifft(Z.prod(axis=0))).real
436#    def _cdf(self, x):
437#    def _ppf(self, q):
438#    def _sf(self, x):
439#    def _isf(self, q):
440#    def _stats(self):
441#    def _entropy(self):
442
443voigt = voigt_gen(name='voigt',shapes='ts,g')
444       
445# Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
446
447class fcjde_gen(st.rv_continuous):
448    """
449    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
450    Ref: J. Appl. Cryst. (1994) 27, 892-900.
451    Parameters
452    -----------------------------------------
453    x: array -1 to 1
454    t: 2-theta position of peak
455    s: sum(S/L,H/L); S: sample height, H: detector opening,
456        L: sample to detector opening distance
457    dx: 2-theta step size in deg
458    Result for fcj.pdf
459    -----------------------------------------
460    T = x*dx+t
461    s = S/L+H/L
462    if x < 0:
463        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
464    if x >= 0:
465        fcj.pdf = 0   
466    """
467    def _pdf(self,x,t,s,dx):
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            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
732            if not GoOn:
733                return -M           #abort!!
734        return M
735       
736    x,y,w,yc,yb,yd = data               #these are numpy arrays!
737    xBeg = np.searchsorted(x,Limits[0])
738    xFin = np.searchsorted(x,Limits[1])
739    bakType,bakDict,bakVary = SetBackgroundParms(Background)
740    dataType,insDict,insVary = SetInstParms(Inst)
741    peakDict,peakVary = setPeaksParms(Peaks)
742    parmDict = {}
743    parmDict.update(bakDict)
744    parmDict.update(insDict)
745    parmDict.update(peakDict)
746    varyList = bakVary+insVary+peakVary
747    while True:
748        begin = time.time()
749        values =  np.array(ValuesOut(parmDict, varyList))
750        if FitPgm == 'LSQ':
751            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
752            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
753            screenSize = wx.ClientDisplayRect()
754            Size = dlg.GetSize()
755            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
756            try:
757                result = so.leastsq(errPeakProfile,values,
758                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg),full_output=True)
759            finally:
760                dlg.Destroy()
761            runtime = time.time()-begin   
762            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
763            print "%s%8.3f%s " % ('fitpeak time =',runtime,'s')
764            ValuesIn(parmDict, varyList, result[0])
765            chisq = np.sum(errPeakProfile(result[0],x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,0)**2)
766            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
767            GOF = chisq/(xFin-xBeg-len(varyList))
768            print "%s%7.2f%s%12.6g%s%6.2f" % ('Rwp = ',Rwp,'%, chi**2 = ',chisq,' reduced chi**2 = ',GOF)
769            try:
770                sig = np.sqrt(np.diag(result[1])*GOF)
771                if np.any(np.isnan(sig)):
772                    print '*** Least squares aborted - some invalid esds possible ***'
773                break                   #refinement succeeded - finish up!
774            except ValueError:          #result[1] is None on singular matrix
775                print 'Refinement failed - singular matrix'
776                Ipvt = result[2]['ipvt']
777                for i,ipvt in enumerate(Ipvt):
778                    if not np.sum(result[2]['fjac'],axis=1)[i]:
779                        print 'Removing parameter: ',varyList[ipvt-1]
780                        del(varyList[ipvt-1])
781                        break
782        elif FitPgm == 'BFGS':
783            print 'Other program here'
784            return
785       
786    sigDict = dict(zip(varyList,sig))
787    yb[xBeg:xFin] = getBackground(parmDict,bakType,x[xBeg:xFin])
788    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
789    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
790    GetBackgroundParms(parmDict,Background)
791    BackgroundPrint(Background,sigDict)
792    GetInstParms(parmDict,Inst,varyList,Peaks)
793    InstPrint(Inst,sigDict)
794    GetPeaksParms(parmDict,Peaks,varyList)   
795    PeaksPrint(parmDict,sigDict,varyList)
796   
797def CalcPDF(data,inst,xydata):
798    auxPlot = []
799    import copy
800    import scipy.fftpack as ft
801    #subtract backgrounds - if any
802    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
803    if data['Sample Bkg.']['Name']:
804        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
805            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
806    if data['Container']['Name']:
807        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
808        if data['Container Bkg.']['Name']:
809            xycontainer += (xydata['Container Bkg.'][1][1]+
810                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
811        xydata['IofQ'][1][1] += xycontainer
812    #get element data & absorption coeff.
813    ElList = data['ElList']
814    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
815    #Apply angle dependent corrections
816    Tth = xydata['Sample'][1][0]
817    dt = (Tth[1]-Tth[0])
818    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
819    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
820    if data['DetType'] == 'Image plate':
821        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
822    XY = xydata['IofQ'][1]   
823    #convert to Q
824    hc = 12.397639
825    if 'Lam' in inst:
826        wave = inst['Lam']
827    else:
828        wave = inst['Lam1']
829    keV = hc/wave
830    minQ = npT2q(Tth[0],wave)
831    maxQ = npT2q(Tth[-1],wave)   
832    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
833    dq = Qpoints[1]-Qpoints[0]
834    XY[0] = npT2q(XY[0],wave)   
835#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
836    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
837    Qdata = T(Qpoints)
838   
839    qLimits = data['QScaleLim']
840    minQ = np.searchsorted(Qpoints,qLimits[0])
841    maxQ = np.searchsorted(Qpoints,qLimits[1])
842    newdata = []
843    xydata['IofQ'][1][0] = Qpoints
844    xydata['IofQ'][1][1] = Qdata
845    for item in xydata['IofQ'][1]:
846        newdata.append(item[:maxQ])
847    xydata['IofQ'][1] = newdata
848   
849
850    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
851    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
852    Q = xydata['SofQ'][1][0]
853    ruland = Ruland(data['Ruland'],wave,Q,CF)
854#    auxPlot.append([Q,ruland,'Ruland'])     
855    CF *= ruland
856#    auxPlot.append([Q,CF,'CF-Corr'])
857    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
858    xydata['SofQ'][1][1] *= scale
859    xydata['SofQ'][1][1] -= CF
860    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
861    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
862    xydata['SofQ'][1][1] *= scale
863   
864    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
865    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
866    if data['Lorch']:
867        xydata['FofQ'][1][1] *= LorchWeight(Q)
868   
869    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
870    nR = len(xydata['GofR'][1][1])
871    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
872    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
873   
874       
875    return auxPlot
876       
877#testing data
878NeedTestData = True
879def TestData():
880#    global NeedTestData
881    NeedTestData = False
882    global bakType
883    bakType = 'chebyschev'
884    global xdata
885    xdata = np.linspace(4.0,40.0,36000)
886    global parmDict0
887    parmDict0 = {
888        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
889        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
890        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
891        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
892        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
893        'Back0':5.384,'Back1':-0.015,'Back2':.004,
894        }
895    global parmDict1
896    parmDict1 = {
897        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
898        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
899        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
900        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
901        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
902        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
903        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
904        'Back0':36.897,'Back1':-0.508,'Back2':.006,
905        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
906        }
907    global varyList
908    varyList = []
909
910def test0():
911    if NeedTestData: TestData()
912    msg = 'test '
913    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
914    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
915    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
916    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
917    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
918    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
919    time0 = time.time()
920    for i in range(100):
921        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
922    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
923   
924
925if __name__ == '__main__':
926    import GSASIItestplot as plot
927    global plotter
928    plotter = plot.PlotNotebook()
929    test0()
930    print "OK"
931    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.