source: trunk/GSASIIpwd.py @ 375

Last change on this file since 375 was 375, checked in by vondreele, 11 years ago

GSASIIgrid.py - add Pawley estimator to menu
GSASIIlattice.py - fix error in CosSinAngle?
GSASIIphsGUI.py - add Pawley estimator & modify size mustrain defaults
GSASIIpwd.py - modify polarization deriv
GSASIIstruct.py - modify refl index issue

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