source: trunk/GSASIIpwd.py @ 353

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

modify psvfcj.for to use sh/l only
add derivative routine to pypowder.for
GSASIIgrid.py - new LS controls for derivative type
GSASIIlattice.py - work on docs as per P Jemian
GSASIIpwd.py - major mods to use FCJpsvoight fortran code & derivatives
GSASIIpwdGUI.py - mod to use LS controls
File size: 43.3 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2011-04-20 13:09:53 -0500 (Wed, 20 Apr 2011) $
6# $Author: vondreele $
7# $Revision: 267 $
8# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 267 2011-04-20 18:09:53Z vondreele $
10########### SVN repository information ###################
11import sys
12import math
13import wx
14import time
15
16import numpy as np
17import scipy as sp
18import numpy.linalg as nl
19from numpy.fft import ifft, fft, fftshift
20import scipy.interpolate as si
21import scipy.stats as st
22import scipy.optimize as so
23
24import GSASIIpath
25import GSASIIplot as G2plt
26import GSASIIlattice as G2lat
27import GSASIIElem as G2elem
28import GSASIIgrid as G2gd
29import GSASIIIO as G2IO
30import pypowder as pyd
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 Transmission(Geometry,Abs,Diam):
96#Calculate sample transmission
97#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
98#   Abs: absorption coeff in cm-1
99#   Diam: sample thickness/diameter in mm
100    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
101        MuR = Abs*Diam/20.0
102        if MuR <= 3.0:
103            T0 = 16/(3.*math.pi)
104            T1 = -0.045780
105            T2 = -0.02489
106            T3 = 0.003045
107            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
108            if T < -20.:
109                return 2.06e-9
110            else:
111                return math.exp(T)
112        else:
113            T1 = 1.433902
114            T2 = 0.013869+0.337894
115            T3 = 1.933433+1.163198
116            T4 = 0.044365-0.04259
117            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
118            return T/100.
119    elif 'plate' in Geometry:
120        MuR = Abs*Diam/10.
121        return math.exp(-MuR)
122    elif 'Bragg' in Geometry:
123        return 0.0
124
125def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
126#Calculate sample absorption
127#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
128#   Abs: absorption coeff in cm-1
129#   Diam: sample thickness/diameter in mm
130#   Tth: 2-theta scattering angle - can be numpy array
131#   Phi: flat plate tilt angle - future
132#   Psi: flat plate tilt axis - future
133    Sth2 = npsind(Tth/2.0)**2
134    Cth2 = 1.-Sth2
135    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
136        MuR = Abs*Diam/20.0
137        if MuR < 3.0:
138            T0 = 16.0/(3*np.pi)
139            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
140                0.109561*np.sqrt(Sth2)-26.04556
141            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
142                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
143            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
144            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
145            return np.exp(Trns)
146        else:
147            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
148                10.02088*Sth2**3-3.36778*Sth2**4
149            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
150                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
151            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
152                0.13576*np.sqrt(Sth2)+1.163198
153            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
154            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
155            return Trns/100.
156    elif 'Bragg' in Geometry:
157        return 1.0
158    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
159        # and only defined for 2theta < 90
160        MuR = Abs*Diam/10.0
161        T1 = np.exp(-MuR)
162        T2 = np.exp(-MuR/npcosd(Tth))
163        Tb = MuR-MuR/npcosd(Tth)
164        return (T2-T1)/Tb
165    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
166        MuR = Abs*Diam/10.0
167        cth = npcosd(Tth/2.0)
168        return np.exp(-MuR/cth)/cth
169       
170def Polarization(Pola,Tth,Azm=0.0):
171#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
172#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
173#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
174#   Tth: 2-theta scattering angle - can be numpy array
175#       which (if either) of these is "right"?
176#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
177#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
178    return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
179        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
180   
181def Oblique(ObCoeff,Tth):
182    if ObCoeff:
183        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
184    else:
185        return 1.0
186               
187def Ruland(RulCoff,wave,Q,Compton):
188    C = 2.9978e8
189    D = 1.5e-3
190    hmc = 0.024262734687
191    sinth2 = (Q*wave/(4.0*np.pi))**2
192    dlam = (wave**2)*Compton*Q/C
193    dlam_c = 2.0*hmc*sinth2-D*wave**2
194    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
195   
196def LorchWeight(Q):
197    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
198           
199def GetAsfMean(ElList,Sthl2):
200#   Calculate various scattering factor terms for PDF calcs
201#   ElList: element dictionary contains scattering factor coefficients, etc.
202#   Sthl2: numpy array of sin theta/lambda squared values
203#   returns: mean(f^2), mean(f)^2, mean(compton)
204    sumNoAtoms = 0.0
205    FF = np.zeros_like(Sthl2)
206    FF2 = np.zeros_like(Sthl2)
207    CF = np.zeros_like(Sthl2)
208    for El in ElList:
209        sumNoAtoms += ElList[El]['FormulaNo']
210    for El in ElList:
211        el = ElList[El]
212        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
213        cf = G2elem.ComptonFac(el,Sthl2)
214        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
215        FF2 += ff2*el['FormulaNo']/sumNoAtoms
216        CF += cf*el['FormulaNo']/sumNoAtoms
217    return FF2,FF**2,CF
218   
219def GetNumDensity(ElList,Vol):
220    sumNoAtoms = 0.0
221    for El in ElList:
222        sumNoAtoms += ElList[El]['FormulaNo']
223    return sumNoAtoms/Vol
224           
225def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
226    NS = 500
227    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
228    Cgama = np.cos(Gama)[:,np.newaxis]
229    Sgama = np.sin(Gama)[:,np.newaxis]
230    CSgama = 1.0/Sgama
231    Delt = Gama[1]-Gama[0]
232    emc = 7.94e-26
233    Navo = 6.023e23
234    Cth = npcosd(Tth/2.0)
235    CTth = npcosd(Tth)
236    Sth = npcosd(Tth/2.0)
237    STth = npsind(Tth)
238    CSth = 1.0/Sth
239    CSTth = 1.0/STth
240    SCth = 1.0/Cth
241    SCTth = 1.0/CTth
242    if 'Bragg' in Geometry:
243        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
244        Y1 = np.pi
245        Y2 = np.pi/2.0
246        Y3 = 3.*np.pi/8. #3pi/4?
247        W = 1.0/(Sth+np.fabs(Sgama))
248        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
249            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
250        Fac0 = Sth**2*Sgama**2
251        X = Fac0+(Fac0+CTth)**2/2
252        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
253        Z = Cgama**4*Cth**4/2.0
254        E = 2.0*(1.0-a)/(b*Cgama/Cth)
255        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
256        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
257        T1 = np.pi/np.sqrt(F1**2-1.0)
258        T2 = np.pi/np.sqrt(F2**2-1.0)
259        Y4 = T1+T2
260        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
261        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
262        Y7 = (T2-T1)/(F1-F2)
263        YT = F2**2*T2-F1**2*T1
264        Y8 = Y1+YT/(F1-F2)
265        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
266        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
267       
268        Q = np.sum(W*M,axis=0)
269        return Q*G*NS/(NS-1.)
270#
271#      cos2th=2.0d*costh^2 - 1.0d
272#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
273#      Y1=3.1415926d
274#      Y2=Y1*0.5d
275#      Y3=Y2*0.75d
276#      for i=1,num_steps-1 do begin
277#         cosgama=double(cos(gama[i]))
278#         singama=double(sin(gama[i]))
279#         cscgama=1.0d / singama
280#
281#         W=1.0d/(sinth+abs(singama))
282#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
283#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
284#
285#         factor0=sinth^2*singama^2
286#         X=factor0+(factor0+cos2th)^2/2.0d
287#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
288#         Z=cosgama^4/2.0d*costh^4
289#         E=2.0d*(1.0-a)/b/cosgama/costh
290#
291#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
292#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
293#
294#         T1=3.14159/sqrt(F1^2-1.0d)
295#         T2=3.14159/sqrt(F2^2-1.0d)
296#         Y4=T1+T2
297#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
298#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
299#         Y7=(T2-T1)/(F1-F2)
300#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
301#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
302#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
303#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
304#
305#         Q=Q+W*M
306#
307#      endfor
308#      Q=double(num_steps)/Double(num_steps-1)*Q*G
309#      end
310    elif 'Cylinder' in Geometry:
311        Q = 0.
312    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
313        EMA = np.exp(-MuT*(1.0-SCTth))
314        Fac1 = (1.-EMA)/(1.0-SCTth)
315        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
316        Fac0 = Cgama/(1-Sgama*SCTth)
317        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
318        Fac0 = Cgama/(1.0+Sgama*SCTth)
319        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
320        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
321        Y = Sgama**3*Cgama*STth*CTth
322        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
323        Fac2 = 1.0/(b*Cgama*STth)
324        U = 2.0*(1.0-a)*Fac2
325        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
326        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)
327        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
328        Y = -Y
329        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)
330        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
331        return Q*G*NS/(NS-1.)
332    elif 'Tilting' in Geometry:
333        EMA = np.exp(-MuT*SCth)
334        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
335#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
336#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
337#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
338#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
339        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
340#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
341        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
342#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
343        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
344        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
345        Z = 0.5*Cgama**4*Sth**4
346#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
347#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
348#
349#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
350#
351#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
352#
353        AlP = 2.0+b*(1.0-Cth*Sgama)
354        AlM = 2.0+b*(1.0+Cth*Sgama)
355#          alphaplus = 2.0 + b*(1.0 - costh*singama)
356#          alphaminus = 2.0 + b*(1.0 + costh*singama)
357        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
358        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
359#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
360#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
361        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)* \
362            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
363            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
364#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
365#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
366#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
367#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
368#                   cosgama*cosgama/b/b*sinth*sinth))
369        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)* \
370            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
371            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
372#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
373#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
374#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
375#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
376#                   cosgama*cosgama/b/b*sinth*sinth))
377        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
378        return Q*G*NS/(NS-1.)
379#       expmutsecth = Exp(-1.0*mut*secth)
380#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
381#       for i=1, num_steps-1 do begin
382#          cosgama=double(cos(gama[i]))
383#          singama=double(sin(gama[i]))
384#          cscgama=1.0d/singama
385#
386#
387#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
388#
389#
390#
391#
392#    ;               print, a,b
393#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
394#          Q=Q+ Wplus*Mplus + Wminus*Mminus
395#      endfor
396#      Q=double(num_steps)/double(num_steps-1)*Q*G
397#   ;   print, min(q), max(q)
398#      end
399
400def MultiScattering(Geometry,ElList,Tth):
401    BN = BD = 0.0
402    Amu = 0.0
403    for El in ElList:
404        el = ElList[El]
405        BN += el['Z']*el['FormulaNo']
406        BD += el['FormulaNo']
407        Amu += el['FormulaNo']*el['mu']
408       
409#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
410
411np.seterr(divide='ignore')
412
413# Normal distribution
414
415# loc = mu, scale = std
416_norm_pdf_C = 1./math.sqrt(2*math.pi)
417class norm_gen(st.rv_continuous):
418       
419    def pdf(self,x,*args,**kwds):
420        loc,scale=kwds['loc'],kwds['scale']
421        x = (x-loc)/scale
422        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
423       
424norm = norm_gen(name='norm',longname='A normal',extradoc="""
425
426Normal distribution
427
428The location (loc) keyword specifies the mean.
429The scale (scale) keyword specifies the standard deviation.
430
431normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
432""")
433
434## Cauchy
435
436# median = loc
437
438class cauchy_gen(st.rv_continuous):
439
440    def pdf(self,x,*args,**kwds):
441        loc,scale=kwds['loc'],kwds['scale']
442        x = (x-loc)/scale
443        return 1.0/np.pi/(1.0+x*x) / scale
444       
445cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
446
447Cauchy distribution
448
449cauchy.pdf(x) = 1/(pi*(1+x**2))
450
451This is the t distribution with one degree of freedom.
452""")
453   
454   
455#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
456
457
458class fcjde_gen(st.rv_continuous):
459    """
460    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
461    Ref: J. Appl. Cryst. (1994) 27, 892-900.
462    Parameters
463    -----------------------------------------
464    x: array -1 to 1
465    t: 2-theta position of peak
466    s: sum(S/L,H/L); S: sample height, H: detector opening,
467        L: sample to detector opening distance
468    dx: 2-theta step size in deg
469    Result for fcj.pdf
470    -----------------------------------------
471    T = x*dx+t
472    s = S/L+H/L
473    if x < 0:
474        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
475    if x >= 0:
476        fcj.pdf = 0   
477    """
478    def _pdf(self,x,t,s,dx):
479        T = dx*x+t
480        ax2 = abs(npcosd(T))
481        ax = ax2**2
482        bx = npcosd(t)**2
483        bx = np.where(ax>bx,bx,ax)
484        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
485        fx = np.where(fx > 0.,fx,0.0)
486        return fx
487             
488    def pdf(self,x,*args,**kwds):
489        loc=kwds['loc']
490        return self._pdf(x-loc,*args)
491       
492fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
493               
494def getWidths(pos,sig,gam,shl):
495    widths = [np.sqrt(sig)/100.,gam/200.]
496    fwhm = 2.355*widths[0]+2.*widths[1]
497    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
498    fmax = 15.0*fwhm
499    if pos > 90:
500        fmin,fmax = [fmax,fmin]         
501    return widths,fmin,fmax
502               
503def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
504    DX = xdata[1]-xdata[0]
505    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
506    x = np.linspace(pos-fmin,pos+fmin,256)
507    dx = x[1]-x[0]
508    Norm = norm.pdf(x,loc=pos,scale=widths[0])
509    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
510    arg = [pos,shl/57.2958,dx,]
511    FCJ = fcjde.pdf(x,*arg,loc=pos)
512    if len(np.nonzero(FCJ)[0])>5:
513        z = np.column_stack([Norm,Cauchy,FCJ]).T
514        Z = fft(z)
515        Df = ifft(Z.prod(axis=0)).real
516    else:
517        z = np.column_stack([Norm,Cauchy]).T
518        Z = fft(z)
519        Df = fftshift(ifft(Z.prod(axis=0))).real
520    Df /= np.sum(Df)
521    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
522    return intens*Df(xdata)*DX/dx
523
524def getBackground(pfx,parmDict,bakType,xdata):
525    yb = np.zeros_like(xdata)
526    if bakType == 'chebyschev':
527        iBak = 0
528        while True:
529            key = pfx+'Back:'+str(iBak)
530            try:
531                yb += parmDict[key]*(xdata-xdata[0])**iBak
532                iBak += 1
533            except KeyError:
534                break
535    return yb
536   
537def getBackgroundDerv(pfx,parmDict,bakType,xdata):
538    dydb = []
539    if bakType == 'chebyschev':
540        iBak = 0
541        while True:
542            if pfx+'Back:'+str(iBak) in parmDict:
543                dydb.append((xdata-xdata[0])**iBak)
544                iBak += 1
545            else:
546                break
547    return dydb
548
549#use old fortran routine
550def getFCJVoigt3(pos,sig,gam,shl,xdata):
551   
552    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
553    Df /= np.sum(Df)
554    return Df
555
556def getdFCJVoigt3(pos,sig,gam,shl,xdata):
557   
558    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) #might have to make these numpy arrays?
559    sumDf = np.sum(Df)
560    return Df/sumDf,dFdp,dFds,dFdg,dFdsh
561   
562
563def getPeakProfile(parmDict,xdata,varyList,bakType):
564   
565    yb = getBackground('',parmDict,bakType,xdata)
566    yc = np.zeros_like(yb)
567    dx = xdata[1]-xdata[0]
568    U = parmDict['U']
569    V = parmDict['V']
570    W = parmDict['W']
571    X = parmDict['X']
572    Y = parmDict['Y']
573    shl = max(parmDict['SH/L'],0.002)
574    Ka2 = False
575    if 'Lam1' in parmDict.keys():
576        Ka2 = True
577        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
578        kRatio = parmDict['I(L2)/I(L1)']
579    iPeak = 0
580    while True:
581        try:
582            pos = parmDict['pos'+str(iPeak)]
583            intens = parmDict['int'+str(iPeak)]
584            sigName = 'sig'+str(iPeak)
585            if sigName in varyList:
586                sig = parmDict[sigName]
587            else:
588                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
589            sig = max(sig,0.001)          #avoid neg sigma
590            gamName = 'gam'+str(iPeak)
591            if gamName in varyList:
592                gam = parmDict[gamName]
593            else:
594                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
595            gam = max(gam,0.001)             #avoid neg gamma
596            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
597            iBeg = np.searchsorted(xdata,pos-fmin)
598            lenX = len(xdata)
599            if not iBeg:
600                iFin = np.searchsorted(xdata,pos+fmin)
601            elif iBeg == lenX:
602                iFin = iBeg
603            else:
604                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
605            if not iBeg+iFin:       #peak below low limit
606                iPeak += 1
607                continue
608            elif not iBeg-iFin:     #peak above high limit
609                return yb+yc
610            yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
611            if Ka2:
612                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
613                kdelt = int((pos2-pos)/dx)               
614                iBeg = min(lenX,iBeg+kdelt)
615                iFin = min(lenX,iFin+kdelt)
616                if iBeg-iFin:
617                    yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
618            iPeak += 1
619        except KeyError:        #no more peaks to process
620            return yb+yc
621           
622def getPeakProfileDerv(parmDict,xdata,varyList,bakType):
623# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
624    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
625    if 'Back:0' in varyList:            #background derivs are in front if present
626        dMdb = getBackgroundDerv('',parmDict,bakType,xdata)
627        dMdv[0:len(dMdb)] = dMdb
628       
629    dx = xdata[1]-xdata[0]
630    U = parmDict['U']
631    V = parmDict['V']
632    W = parmDict['W']
633    X = parmDict['X']
634    Y = parmDict['Y']
635    shl = max(parmDict['SH/L'],0.002)
636    Ka2 = False
637    if 'Lam1' in parmDict.keys():
638        Ka2 = True
639        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
640        kRatio = parmDict['I(L2)/I(L1)']
641    iPeak = 0
642    while True:
643        try:
644            pos = parmDict['pos'+str(iPeak)]
645            intens = parmDict['int'+str(iPeak)]
646            sigName = 'sig'+str(iPeak)
647            tanth = tand(pos/2.0)
648            costh = cosd(pos/2.0)
649            if sigName in varyList:
650                sig = parmDict[sigName]
651            else:
652                sig = U*tanth**2+V*tanth+W
653                dsdU = tanth**2
654                dsdV = tanth
655                dsdW = 1.0
656            sig = max(sig,0.001)          #avoid neg sigma
657            gamName = 'gam'+str(iPeak)
658            if gamName in varyList:
659                gam = parmDict[gamName]
660            else:
661                gam = X/costh+Y*tanth
662                dgdX = 1.0/costh
663                dgdY = tanth
664            gam = max(gam,0.001)             #avoid neg gamma
665            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
666            iBeg = np.searchsorted(xdata,pos-fmin)
667            lenX = len(xdata)
668            if not iBeg:
669                iFin = np.searchsorted(xdata,pos+fmin)
670            elif iBeg == lenX:
671                iFin = iBeg
672            else:
673                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
674            if not iBeg+iFin:       #peak below low limit
675                iPeak += 1
676                continue
677            elif not iBeg-iFin:     #peak above high limit
678                break
679            dMdpk = np.zeros(shape=(6,len(xdata)))
680            dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
681            for i in range(5):
682                dMdpk[i][iBeg:iFin] = 100.*dx*intens*dMdipk[i]
683            dMdpk[0][iBeg:iFin] = dMdipk[0]
684            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
685            if Ka2:
686                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
687                kdelt = int((pos2-pos)/dx)               
688                iBeg = min(lenX,iBeg+kdelt)
689                iFin = min(lenX,iFin+kdelt)
690                if iBeg-iFin:
691                    dMdipk = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
692                    for i in range(5):
693                        dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk[i]
694                    dMdpk[0][iBeg:iFin] += kRatio*dMdipk[0]
695                    dMdpk[5][iBeg:iFin] += dMdipk[0]
696                    dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
697            for parmName in ['pos','int','sig','gam']:
698                try:
699                    idx = varyList.index(parmName+str(iPeak))
700                    dMdv[idx] = dervDict[parmName]
701                except ValueError:
702                    pass
703            if 'U' in varyList:
704                dMdv[varyList.index('U')] += dsdU*dervDict['sig']
705            if 'V' in varyList:
706                dMdv[varyList.index('V')] += dsdV*dervDict['sig']
707            if 'W' in varyList:
708                dMdv[varyList.index('W')] += dsdW*dervDict['sig']
709            if 'X' in varyList:
710                dMdv[varyList.index('X')] += dgdX*dervDict['gam']
711            if 'Y' in varyList:
712                dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
713            if 'SH/L' in varyList:
714                dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
715            if 'I(L2)/I(L1)' in varyList:
716                dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
717            iPeak += 1
718        except KeyError:        #no more peaks to process
719            break
720    return dMdv
721       
722def Dict2Values(parmdict, varylist):
723    '''Use before call to leastsq to setup list of values for the parameters
724    in parmdict, as selected by key in varylist'''
725    return [parmdict[key] for key in varylist] 
726   
727def Values2Dict(parmdict, varylist, values):
728    ''' Use after call to leastsq to update the parameter dictionary with
729    values corresponding to keys in varylist'''
730    parmdict.update(zip(varylist,values))
731   
732def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None):
733   
734    def SetBackgroundParms(Background):
735        bakType,bakFlag = Background[:2]
736        backVals = Background[3:]
737        backNames = ['Back:'+str(i) for i in range(len(backVals))]
738        if bakFlag: #returns backNames as varyList = backNames
739            return bakType,dict(zip(backNames,backVals)),backNames
740        else:       #no background varied; varyList = []
741            return bakType,dict(zip(backNames,backVals)),[]
742       
743    def GetBackgroundParms(parmList,Background):
744        iBak = 0
745        while True:
746            try:
747                bakName = 'Back:'+str(iBak)
748                Background[iBak+3] = parmList[bakName]
749                iBak += 1
750            except KeyError:
751                break
752               
753    def BackgroundPrint(Background,sigDict):
754        if Background[1]:
755            print 'Background coefficients for',Background[0],'function'
756            ptfmt = "%12.5f"
757            ptstr =  'values:'
758            sigstr = 'esds  :'
759            for i,back in enumerate(Background[3:]):
760                ptstr += ptfmt % (back)
761                sigstr += ptfmt % (sigDict['Back:'+str(i)])
762            print ptstr
763            print sigstr
764        else:
765            print 'Background not refined'
766           
767    def SetInstParms(Inst):
768        insVals,insFlags,insNames = Inst[1:4]
769        dataType = insVals[0]
770        insVary = []
771        for i,flag in enumerate(insFlags):
772            if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
773                insVary.append(insNames[i])
774        instDict = dict(zip(insNames,insVals))
775        instDict['X'] = max(instDict['X'],0.01)
776        instDict['Y'] = max(instDict['Y'],0.01)
777        instDict['SH/L'] = max(instDict['SH/L'],0.002)
778        return dataType,instDict,insVary
779       
780    def GetInstParms(parmDict,Inst,varyList,Peaks):
781        instNames = Inst[3]
782        for i,name in enumerate(instNames):
783            Inst[1][i] = parmDict[name]
784        iPeak = 0
785        while True:
786            try:
787                sigName = 'sig'+str(iPeak)
788                pos = parmDict['pos'+str(iPeak)]
789                if sigName not in varyList:
790                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
791                gamName = 'gam'+str(iPeak)
792                if gamName not in varyList:
793                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
794                iPeak += 1
795            except KeyError:
796                break
797       
798    def InstPrint(Inst,sigDict):
799        print 'Instrument Parameters:'
800        ptfmt = "%12.6f"
801        ptlbls = 'names :'
802        ptstr =  'values:'
803        sigstr = 'esds  :'
804        instNames = Inst[3][1:]
805        for i,name in enumerate(instNames):
806            ptlbls += "%s" % (name.center(12))
807            ptstr += ptfmt % (Inst[1][i+1])
808            if name in sigDict:
809                sigstr += ptfmt % (sigDict[name])
810            else:
811                sigstr += 12*' '
812        print ptlbls
813        print ptstr
814        print sigstr
815
816    def SetPeaksParms(Peaks):
817        peakNames = []
818        peakVary = []
819        peakVals = []
820        names = ['pos','int','sig','gam']
821        for i,peak in enumerate(Peaks):
822            for j in range(4):
823                peakVals.append(peak[2*j])
824                parName = names[j]+str(i)
825                peakNames.append(parName)
826                if peak[2*j+1]:
827                    peakVary.append(parName)
828        return dict(zip(peakNames,peakVals)),peakVary
829               
830    def GetPeaksParms(parmDict,Peaks,varyList):
831        names = ['pos','int','sig','gam']
832        for i,peak in enumerate(Peaks):
833            for j in range(4):
834                pos = parmDict['pos'+str(i)]
835                parName = names[j]+str(i)
836                if parName in varyList:
837                    peak[2*j] = parmDict[parName]
838                elif 'sig' in parName:
839                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
840                elif 'gam' in parName:
841                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
842                       
843    def PeaksPrint(parmDict,sigDict,varyList):
844        print 'Peak coefficients:'
845        names = ['pos','int','sig','gam']
846        head = 15*' '
847        for name in names:
848            head += name.center(12)+'esd'.center(12)
849        print head
850        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
851        for i,peak in enumerate(Peaks):
852            ptstr =  ':'
853            for j in range(4):
854                name = names[j]
855                parName = name+str(i)
856                ptstr += ptfmt[name] % (parmDict[parName])
857                if parName in varyList:
858#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
859                    ptstr += ptfmt[name] % (sigDict[parName])
860                else:
861#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
862                    ptstr += 12*' '
863            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
864               
865    def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):
866        parmdict.update(zip(varylist,values))
867        return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType)
868           
869    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
870        parmdict.update(zip(varylist,values))
871        M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata)
872        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
873        if dlg:
874            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
875            if not GoOn:
876                return -M           #abort!!
877        return M
878       
879    if controls:
880        Ftol = controls['min dM/M']
881        derivType = controls['deriv type']
882    else:
883        Ftol = 0.0001
884        derivType = 'analytic'
885    if oneCycle:
886        Ftol = 1.0
887    x,y,w,yc,yb,yd = data               #these are numpy arrays!
888    xBeg = np.searchsorted(x,Limits[0])
889    xFin = np.searchsorted(x,Limits[1])
890    bakType,bakDict,bakVary = SetBackgroundParms(Background)
891    dataType,insDict,insVary = SetInstParms(Inst)
892    peakDict,peakVary = SetPeaksParms(Peaks)
893    parmDict = {}
894    parmDict.update(bakDict)
895    parmDict.update(insDict)
896    parmDict.update(peakDict)
897    varyList = bakVary+insVary+peakVary
898    while True:
899        begin = time.time()
900        values =  np.array(Dict2Values(parmDict, varyList))
901        if FitPgm == 'LSQ':
902            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
903            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
904            screenSize = wx.ClientDisplayRect()
905            Size = dlg.GetSize()
906            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
907            try:
908                if derivType == 'analytic':
909                    result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
910                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
911                    ncyc = int(result[2]['nfev']/2)
912                else:
913                    result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
914                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
915                    ncyc = int(result[2]['nfev']/len(varyList))
916            finally:
917                dlg.Destroy()
918            runtime = time.time()-begin   
919            chisq = np.sum(result[2]['fvec']**2)
920            Values2Dict(parmDict, varyList, result[0])
921            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
922            GOF = chisq/(xFin-xBeg-len(varyList))
923            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
924            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
925            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
926            try:
927                sig = np.sqrt(np.diag(result[1])*GOF)
928                if np.any(np.isnan(sig)):
929                    print '*** Least squares aborted - some invalid esds possible ***'
930                break                   #refinement succeeded - finish up!
931            except ValueError:          #result[1] is None on singular matrix
932                print '**** Refinement failed - singular matrix ****'
933                Ipvt = result[2]['ipvt']
934                for i,ipvt in enumerate(Ipvt):
935                    if not np.sum(result[2]['fjac'],axis=1)[i]:
936                        print 'Removing parameter: ',varyList[ipvt-1]
937                        del(varyList[ipvt-1])
938                        break
939        elif FitPgm == 'BFGS':
940            print 'Other program here'
941            return
942       
943    sigDict = dict(zip(varyList,sig))
944    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
945    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
946    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
947    GetBackgroundParms(parmDict,Background)
948    BackgroundPrint(Background,sigDict)
949    GetInstParms(parmDict,Inst,varyList,Peaks)
950    InstPrint(Inst,sigDict)
951    GetPeaksParms(parmDict,Peaks,varyList)   
952    PeaksPrint(parmDict,sigDict,varyList)
953   
954def CalcPDF(data,inst,xydata):
955    auxPlot = []
956    import copy
957    import scipy.fftpack as ft
958    #subtract backgrounds - if any
959    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
960    if data['Sample Bkg.']['Name']:
961        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
962            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
963    if data['Container']['Name']:
964        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
965        if data['Container Bkg.']['Name']:
966            xycontainer += (xydata['Container Bkg.'][1][1]+
967                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
968        xydata['IofQ'][1][1] += xycontainer
969    #get element data & absorption coeff.
970    ElList = data['ElList']
971    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
972    #Apply angle dependent corrections
973    Tth = xydata['Sample'][1][0]
974    dt = (Tth[1]-Tth[0])
975    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
976    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
977    if data['DetType'] == 'Image plate':
978        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
979    XY = xydata['IofQ'][1]   
980    #convert to Q
981    hc = 12.397639
982    if 'Lam' in inst:
983        wave = inst['Lam']
984    else:
985        wave = inst['Lam1']
986    keV = hc/wave
987    minQ = npT2q(Tth[0],wave)
988    maxQ = npT2q(Tth[-1],wave)   
989    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
990    dq = Qpoints[1]-Qpoints[0]
991    XY[0] = npT2q(XY[0],wave)   
992#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
993    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
994    Qdata = T(Qpoints)
995   
996    qLimits = data['QScaleLim']
997    minQ = np.searchsorted(Qpoints,qLimits[0])
998    maxQ = np.searchsorted(Qpoints,qLimits[1])
999    newdata = []
1000    xydata['IofQ'][1][0] = Qpoints
1001    xydata['IofQ'][1][1] = Qdata
1002    for item in xydata['IofQ'][1]:
1003        newdata.append(item[:maxQ])
1004    xydata['IofQ'][1] = newdata
1005   
1006
1007    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
1008    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
1009    Q = xydata['SofQ'][1][0]
1010    ruland = Ruland(data['Ruland'],wave,Q,CF)
1011#    auxPlot.append([Q,ruland,'Ruland'])     
1012    CF *= ruland
1013#    auxPlot.append([Q,CF,'CF-Corr'])
1014    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
1015    xydata['SofQ'][1][1] *= scale
1016    xydata['SofQ'][1][1] -= CF
1017    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
1018    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
1019    xydata['SofQ'][1][1] *= scale
1020   
1021    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
1022    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
1023    if data['Lorch']:
1024        xydata['FofQ'][1][1] *= LorchWeight(Q)
1025   
1026    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
1027    nR = len(xydata['GofR'][1][1])
1028    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
1029    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
1030   
1031       
1032    return auxPlot
1033       
1034#testing data
1035NeedTestData = True
1036def TestData():
1037#    global NeedTestData
1038    NeedTestData = False
1039    global bakType
1040    bakType = 'chebyschev'
1041    global xdata
1042    xdata = np.linspace(4.0,40.0,36000)
1043    global parmDict0
1044    parmDict0 = {
1045        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1046        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1047        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1048        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1049        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1050        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1051        }
1052    global parmDict1
1053    parmDict1 = {
1054        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1055        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1056        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1057        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1058        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1059        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1060        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1061        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1062        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1063        }
1064    global parmDict2
1065    parmDict2 = {
1066        'pos0':5.7,'int0':10.0,'sig0':0.5,'gam0':1.0,
1067        'U':1.,'V':-1.,'W':0.1,'X':0.0,'Y':2.,'SH/L':0.004,
1068        'Back0':5.,'Back1':-0.02,'Back2':.004,
1069        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1070        }
1071    global varyList
1072    varyList = []
1073
1074def test0():
1075    if NeedTestData: TestData()
1076    msg = 'test '
1077    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1078    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1079    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1080    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1081    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1082    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1083   
1084def test1():
1085    if NeedTestData: TestData()
1086    time0 = time.time()
1087    for i in range(100):
1088        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1089    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1090   
1091def test2(name,delt):
1092    if NeedTestData: TestData()
1093    varyList = [name,]
1094    xdata = np.linspace(5.6,5.8,800)
1095    hplot = plotter.add('derivatives test for '+name).gca()
1096    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1097    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1098    parmDict2[name] += delt
1099    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1100    hplot.plot(xdata,(y1-y0)/delt,'r+')
1101   
1102   
1103
1104if __name__ == '__main__':
1105    import GSASIItestplot as plot
1106    global plotter
1107    plotter = plot.PlotNotebook()
1108#    test0()
1109    test2('I(L2)/I(L1)',.0001)
1110    print "OK"
1111    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.