source: trunk/GSASIIpwd.py @ 335

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

Add delete Pawley list; scratch "leBail" from things
mods to texture display stuff
refactor instrument parameters GUI as sizer based - no grid table
some mods to peak fit output
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            chisq = np.sum(result[2]['fvec']**2)
763            ncyc = int(result[2]['nfev']/len(varyList))
764            ValuesIn(parmDict, varyList, result[0])
765            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
766            GOF = chisq/(xFin-xBeg-len(varyList))
767            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
768            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
769            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
770            try:
771                sig = np.sqrt(np.diag(result[1])*GOF)
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
879NeedTestData = True
880def TestData():
881#    global NeedTestData
882    NeedTestData = False
883    global bakType
884    bakType = 'chebyschev'
885    global xdata
886    xdata = np.linspace(4.0,40.0,36000)
887    global parmDict0
888    parmDict0 = {
889        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
890        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
891        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
892        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
893        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
894        'Back0':5.384,'Back1':-0.015,'Back2':.004,
895        }
896    global parmDict1
897    parmDict1 = {
898        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
899        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
900        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
901        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
902        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
903        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
904        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
905        'Back0':36.897,'Back1':-0.508,'Back2':.006,
906        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
907        }
908    global varyList
909    varyList = []
910
911def test0():
912    if NeedTestData: TestData()
913    msg = 'test '
914    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
915    gplot.plot(xdata,getBackground(parmDict0,bakType,xdata))   
916    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
917    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
918    fplot.plot(xdata,getBackground(parmDict1,bakType,xdata))   
919    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
920    time0 = time.time()
921    for i in range(100):
922        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
923    print '100+6*Ka1-2 peaks=3600 peaks',time.time()-time0
924   
925
926if __name__ == '__main__':
927    import GSASIItestplot as plot
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.