source: trunk/GSASIIpwd.py @ 797

Last change on this file since 797 was 797, checked in by vondreele, 9 years ago

some fixes to TOF stuff
some fixes to SXC stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 64.3 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2012-11-08 22:05:35 +0000 (Thu, 08 Nov 2012) $
6# $Author: vondreele $
7# $Revision: 797 $
8# $URL: trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 797 2012-11-08 22:05:35Z vondreele $
10########### SVN repository information ###################
11import sys
12import math
13import time
14
15import numpy as np
16import scipy as sp
17import numpy.linalg as nl
18from numpy.fft import ifft, fft, fftshift
19import scipy.interpolate as si
20import scipy.stats as st
21import scipy.optimize as so
22
23import GSASIIpath
24GSASIIpath.SetVersionNumber("$Revision: 797 $")
25import GSASIIlattice as G2lat
26import GSASIIspc as G2spc
27import GSASIIElem as G2elem
28import GSASIIgrid as G2gd
29import GSASIIIO as G2IO
30import GSASIImath as G2mth
31import pypowder as pyd
32
33# trig functions in degrees
34sind = lambda x: math.sin(x*math.pi/180.)
35asind = lambda x: 180.*math.asin(x)/math.pi
36tand = lambda x: math.tan(x*math.pi/180.)
37atand = lambda x: 180.*math.atan(x)/math.pi
38atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
39cosd = lambda x: math.cos(x*math.pi/180.)
40acosd = lambda x: 180.*math.acos(x)/math.pi
41rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
42#numpy versions
43npsind = lambda x: np.sin(x*np.pi/180.)
44npasind = lambda x: 180.*np.arcsin(x)/math.pi
45npcosd = lambda x: np.cos(x*math.pi/180.)
46npacosd = lambda x: 180.*np.arccos(x)/math.pi
47nptand = lambda x: np.tan(x*math.pi/180.)
48npatand = lambda x: 180.*np.arctan(x)/np.pi
49npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
51npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
52   
53#GSASII pdf calculation routines
54       
55def Transmission(Geometry,Abs,Diam):
56#Calculate sample transmission
57#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
58#   Abs: absorption coeff in cm-1
59#   Diam: sample thickness/diameter in mm
60    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
61        MuR = Abs*Diam/20.0
62        if MuR <= 3.0:
63            T0 = 16/(3.*math.pi)
64            T1 = -0.045780
65            T2 = -0.02489
66            T3 = 0.003045
67            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
68            if T < -20.:
69                return 2.06e-9
70            else:
71                return math.exp(T)
72        else:
73            T1 = 1.433902
74            T2 = 0.013869+0.337894
75            T3 = 1.933433+1.163198
76            T4 = 0.044365-0.04259
77            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
78            return T/100.
79    elif 'plate' in Geometry:
80        MuR = Abs*Diam/10.
81        return math.exp(-MuR)
82    elif 'Bragg' in Geometry:
83        return 0.0
84
85def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
86#Calculate sample absorption
87#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
88#   MuR: absorption coeff * sample thickness/2 or radius
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        if MuR < 3.0:
96            T0 = 16.0/(3*np.pi)
97            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
98                0.109561*np.sqrt(Sth2)-26.04556
99            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
100                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
101            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
102            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
103            return np.exp(Trns)
104        else:
105            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
106                10.02088*Sth2**3-3.36778*Sth2**4
107            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
108                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
109            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
110                0.13576*np.sqrt(Sth2)+1.163198
111            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
112            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
113            return Trns/100.
114    elif 'Bragg' in Geometry:
115        return 1.0
116    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
117        # and only defined for 2theta < 90
118        MuT = 2.*MuR
119        T1 = np.exp(-MuT)
120        T2 = np.exp(-MuT/npcosd(Tth))
121        Tb = MuT-MuT/npcosd(Tth)
122        return (T2-T1)/Tb
123    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
124        MuT = 2.*MuR
125        cth = npcosd(Tth/2.0)
126        return np.exp(-MuT/cth)/cth
127       
128def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
129    dA = 0.001
130    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
131    if MuR:
132        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
133        return (AbsP-AbsM)/(2.0*dA)
134    else:
135        return (AbsP-1.)/dA
136       
137def Polarization(Pola,Tth,Azm=0.0):
138    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
139    input:
140        Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
141        Azm: azimuthal angle e.g. 0.0 in plane of polarization
142        Tth: 2-theta scattering angle - can be numpy array
143            which (if either) of these is "right"?
144    return:
145        pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
146            (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
147        dpdPola: derivative needed for least squares
148    """
149    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
150        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
151    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
152    return pola,dpdPola
153   
154def Oblique(ObCoeff,Tth):
155    if ObCoeff:
156        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
157    else:
158        return 1.0
159               
160def Ruland(RulCoff,wave,Q,Compton):
161    C = 2.9978e8
162    D = 1.5e-3
163    hmc = 0.024262734687
164    sinth2 = (Q*wave/(4.0*np.pi))**2
165    dlam = (wave**2)*Compton*Q/C
166    dlam_c = 2.0*hmc*sinth2-D*wave**2
167    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
168   
169def LorchWeight(Q):
170    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
171           
172def GetAsfMean(ElList,Sthl2):
173#   Calculate various scattering factor terms for PDF calcs
174#   ElList: element dictionary contains scattering factor coefficients, etc.
175#   Sthl2: numpy array of sin theta/lambda squared values
176#   returns: mean(f^2), mean(f)^2, mean(compton)
177    sumNoAtoms = 0.0
178    FF = np.zeros_like(Sthl2)
179    FF2 = np.zeros_like(Sthl2)
180    CF = np.zeros_like(Sthl2)
181    for El in ElList:
182        sumNoAtoms += ElList[El]['FormulaNo']
183    for El in ElList:
184        el = ElList[El]
185        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
186        cf = G2elem.ComptonFac(el,Sthl2)
187        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
188        FF2 += ff2*el['FormulaNo']/sumNoAtoms
189        CF += cf*el['FormulaNo']/sumNoAtoms
190    return FF2,FF**2,CF
191   
192def GetNumDensity(ElList,Vol):
193    sumNoAtoms = 0.0
194    for El in ElList:
195        sumNoAtoms += ElList[El]['FormulaNo']
196    return sumNoAtoms/Vol
197           
198#def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
199#    NS = 500
200#    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
201#    Cgama = np.cos(Gama)[:,np.newaxis]
202#    Sgama = np.sin(Gama)[:,np.newaxis]
203#    CSgama = 1.0/Sgama
204#    Delt = Gama[1]-Gama[0]
205#    emc = 7.94e-26
206#    Navo = 6.023e23
207#    Cth = npcosd(Tth/2.0)
208#    CTth = npcosd(Tth)
209#    Sth = npcosd(Tth/2.0)
210#    STth = npsind(Tth)
211#    CSth = 1.0/Sth
212#    CSTth = 1.0/STth
213#    SCth = 1.0/Cth
214#    SCTth = 1.0/CTth
215#    if 'Bragg' in Geometry:
216#        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
217#        Y1 = np.pi
218#        Y2 = np.pi/2.0
219#        Y3 = 3.*np.pi/8. #3pi/4?
220#        W = 1.0/(Sth+np.fabs(Sgama))
221#        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
222#            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
223#        Fac0 = Sth**2*Sgama**2
224#        X = Fac0+(Fac0+CTth)**2/2
225#        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
226#        Z = Cgama**4*Cth**4/2.0
227#        E = 2.0*(1.0-a)/(b*Cgama/Cth)
228#        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
229#        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
230#        T1 = np.pi/np.sqrt(F1**2-1.0)
231#        T2 = np.pi/np.sqrt(F2**2-1.0)
232#        Y4 = T1+T2
233#        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
234#        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
235#        Y7 = (T2-T1)/(F1-F2)
236#        YT = F2**2*T2-F1**2*T1
237#        Y8 = Y1+YT/(F1-F2)
238#        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
239#        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
240#       
241#        Q = np.sum(W*M,axis=0)
242#        return Q*G*NS/(NS-1.)
243##
244##      cos2th=2.0d*costh^2 - 1.0d
245##      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
246##      Y1=3.1415926d
247##      Y2=Y1*0.5d
248##      Y3=Y2*0.75d
249##      for i=1,num_steps-1 do begin
250##         cosgama=double(cos(gama[i]))
251##         singama=double(sin(gama[i]))
252##         cscgama=1.0d / singama
253##
254##         W=1.0d/(sinth+abs(singama))
255##         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
256##                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
257##
258##         factor0=sinth^2*singama^2
259##         X=factor0+(factor0+cos2th)^2/2.0d
260##         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
261##         Z=cosgama^4/2.0d*costh^4
262##         E=2.0d*(1.0-a)/b/cosgama/costh
263##
264##         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
265##         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
266##
267##         T1=3.14159/sqrt(F1^2-1.0d)
268##         T2=3.14159/sqrt(F2^2-1.0d)
269##         Y4=T1+T2
270##         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
271##         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
272##         Y7=(T2-T1)/(F1-F2)
273##         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
274##         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
275##         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
276##                      (X*Y7+Y*Y8+Z*Y9))*cosgama
277##
278##         Q=Q+W*M
279##
280##      endfor
281##      Q=double(num_steps)/Double(num_steps-1)*Q*G
282##      end
283#    elif 'Cylinder' in Geometry:
284#        Q = 0.
285#    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
286#        EMA = np.exp(-MuT*(1.0-SCTth))
287#        Fac1 = (1.-EMA)/(1.0-SCTth)
288#        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
289#        Fac0 = Cgama/(1-Sgama*SCTth)
290#        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
291#        Fac0 = Cgama/(1.0+Sgama*SCTth)
292#        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
293#        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
294#        Y = Sgama**3*Cgama*STth*CTth
295#        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
296#        Fac2 = 1.0/(b*Cgama*STth)
297#        U = 2.0*(1.0-a)*Fac2
298#        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
299#        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)
300#        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
301#        Y = -Y
302#        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)
303#        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
304#        return Q*G*NS/(NS-1.)
305#    elif 'Tilting' in Geometry:
306#        EMA = np.exp(-MuT*SCth)
307#        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
308##          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
309##                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
310##          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
311##                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
312#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
313##        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
314#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
315##        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
316#        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
317#        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
318#        Z = 0.5*Cgama**4*Sth**4
319##          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
320##                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
321##
322##          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
323##
324##          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
325##
326#        AlP = 2.0+b*(1.0-Cth*Sgama)
327#        AlM = 2.0+b*(1.0+Cth*Sgama)
328##          alphaplus = 2.0 + b*(1.0 - costh*singama)
329##          alphaminus = 2.0 + b*(1.0 + costh*singama)
330#        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
331#        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
332##          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
333##          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
334#        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)* \
335#            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
336#            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
337##          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
338##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
339##                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
340##                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
341##                   cosgama*cosgama/b/b*sinth*sinth))
342#        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)* \
343#            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
344#            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
345##          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
346##                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
347##                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
348##                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
349##                   cosgama*cosgama/b/b*sinth*sinth))
350#        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
351#        return Q*G*NS/(NS-1.)
352##       expmutsecth = Exp(-1.0*mut*secth)
353##       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
354##       for i=1, num_steps-1 do begin
355##          cosgama=double(cos(gama[i]))
356##          singama=double(sin(gama[i]))
357##          cscgama=1.0d/singama
358##
359##
360##     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
361##
362##
363##
364##
365##    ;               print, a,b
366##  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
367##          Q=Q+ Wplus*Mplus + Wminus*Mminus
368##      endfor
369##      Q=double(num_steps)/double(num_steps-1)*Q*G
370##   ;   print, min(q), max(q)
371##      end
372
373#def MultiScattering(Geometry,ElList,Tth):
374#    BN = BD = 0.0
375#    Amu = 0.0
376#    for El in ElList:
377#        el = ElList[El]
378#        BN += el['Z']*el['FormulaNo']
379#        BD += el['FormulaNo']
380#        Amu += el['FormulaNo']*el['mu']
381       
382def CalcPDF(data,inst,xydata):
383    auxPlot = []
384    import copy
385    import scipy.fftpack as ft
386    #subtract backgrounds - if any
387    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
388    if data['Sample Bkg.']['Name']:
389        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
390            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
391    if data['Container']['Name']:
392        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
393        if data['Container Bkg.']['Name']:
394            xycontainer += (xydata['Container Bkg.'][1][1]+
395                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
396        xydata['IofQ'][1][1] += xycontainer
397    #get element data & absorption coeff.
398    ElList = data['ElList']
399    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
400    #Apply angle dependent corrections
401    Tth = xydata['Sample'][1][0]
402    dt = (Tth[1]-Tth[0])
403    MuR = Abs*data['Diam']/20.0
404    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
405    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
406    if data['DetType'] == 'Image plate':
407        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
408    XY = xydata['IofQ'][1]   
409    #convert to Q
410    hc = 12.397639
411    wave = G2mth.getWave(inst)
412    keV = hc/wave
413    minQ = npT2q(Tth[0],wave)
414    maxQ = npT2q(Tth[-1],wave)   
415    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
416    dq = Qpoints[1]-Qpoints[0]
417    XY[0] = npT2q(XY[0],wave)   
418#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
419    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
420    Qdata = T(Qpoints)
421   
422    qLimits = data['QScaleLim']
423    minQ = np.searchsorted(Qpoints,qLimits[0])
424    maxQ = np.searchsorted(Qpoints,qLimits[1])
425    newdata = []
426    xydata['IofQ'][1][0] = Qpoints
427    xydata['IofQ'][1][1] = Qdata
428    for item in xydata['IofQ'][1]:
429        newdata.append(item[:maxQ])
430    xydata['IofQ'][1] = newdata
431   
432
433    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
434    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
435    Q = xydata['SofQ'][1][0]
436    ruland = Ruland(data['Ruland'],wave,Q,CF)
437#    auxPlot.append([Q,ruland,'Ruland'])     
438    CF *= ruland
439#    auxPlot.append([Q,CF,'CF-Corr'])
440    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
441    xydata['SofQ'][1][1] *= scale
442    xydata['SofQ'][1][1] -= CF
443    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
444    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
445    xydata['SofQ'][1][1] *= scale
446   
447    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
448    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
449    if data['Lorch']:
450        xydata['FofQ'][1][1] *= LorchWeight(Q)
451   
452    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
453    nR = len(xydata['GofR'][1][1])
454    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
455    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
456   
457       
458    return auxPlot
459       
460#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
461
462def factorize(num):
463    ''' Provide prime number factors for integer num
464    Returns dictionary of prime factors (keys) & power for each (data)
465    '''
466    factors = {}
467    orig = num
468
469    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
470    i, sqi = 2, 4
471    while sqi <= num:
472        while not num%i:
473            num /= i
474            factors[i] = factors.get(i, 0) + 1
475
476        sqi += 2*i + 1
477        i += 1
478
479    if num != 1 and num != orig:
480        factors[num] = factors.get(num, 0) + 1
481
482    if factors:
483        return factors
484    else:
485        return {num:1}          #a prime number!
486           
487def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
488    ''' Provide list of optimal data sizes for FFT calculations
489    Input:
490        nmin: minimum data size >= 1
491        nmax: maximum data size > nmin
492        thresh: maximum prime factor allowed
493    Returns:
494        list of data sizes where the maximum prime factor is < thresh
495    ''' 
496    plist = []
497    nmin = max(1,nmin)
498    nmax = max(nmin+1,nmax)
499    for p in range(nmin,nmax):
500        if max(factorize(p).keys()) < thresh:
501            plist.append(p)
502    return plist
503
504np.seterr(divide='ignore')
505
506# Normal distribution
507
508# loc = mu, scale = std
509_norm_pdf_C = 1./math.sqrt(2*math.pi)
510class norm_gen(st.rv_continuous):
511       
512    def pdf(self,x,*args,**kwds):
513        loc,scale=kwds['loc'],kwds['scale']
514        x = (x-loc)/scale
515        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
516       
517norm = norm_gen(name='norm',longname='A normal',extradoc="""
518
519Normal distribution
520
521The location (loc) keyword specifies the mean.
522The scale (scale) keyword specifies the standard deviation.
523
524normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
525""")
526
527## Cauchy
528
529# median = loc
530
531class cauchy_gen(st.rv_continuous):
532
533    def pdf(self,x,*args,**kwds):
534        loc,scale=kwds['loc'],kwds['scale']
535        x = (x-loc)/scale
536        return 1.0/np.pi/(1.0+x*x) / scale
537       
538cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
539
540Cauchy distribution
541
542cauchy.pdf(x) = 1/(pi*(1+x**2))
543
544This is the t distribution with one degree of freedom.
545""")
546   
547   
548#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
549
550
551class fcjde_gen(st.rv_continuous):
552    """
553    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
554    Ref: J. Appl. Cryst. (1994) 27, 892-900.
555    Parameters
556    -----------------------------------------
557    x: array -1 to 1
558    t: 2-theta position of peak
559    s: sum(S/L,H/L); S: sample height, H: detector opening,
560        L: sample to detector opening distance
561    dx: 2-theta step size in deg
562    Result for fcj.pdf
563    -----------------------------------------
564    T = x*dx+t
565    s = S/L+H/L
566    if x < 0:
567        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
568    if x >= 0:
569        fcj.pdf = 0   
570    """
571    def _pdf(self,x,t,s,dx):
572        T = dx*x+t
573        ax2 = abs(npcosd(T))
574        ax = ax2**2
575        bx = npcosd(t)**2
576        bx = np.where(ax>bx,bx,ax)
577        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
578        fx = np.where(fx > 0.,fx,0.0)
579        return fx
580             
581    def pdf(self,x,*args,**kwds):
582        loc=kwds['loc']
583        return self._pdf(x-loc,*args)
584       
585fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
586               
587def getWidthsCW(pos,sig,gam,shl):
588    widths = [np.sqrt(sig)/100.,gam/200.]
589    fwhm = 2.355*widths[0]+2.*widths[1]
590    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
591    fmax = 15.0*fwhm
592    if pos > 90:
593        fmin,fmax = [fmax,fmin]         
594    return widths,fmin,fmax
595   
596def getWidthsTOF(pos,alp,bet,sig,gam):
597    lnf = 3.3      # =log(0.001/2)
598    widths = [np.sqrt(sig),gam]
599    fwhm = 2.355*widths[0]+2.*widths[1]
600    fmin = 10.*fwhm*(1.+1./alp)   
601    fmax = 10.*fwhm*(1.+1./bet)
602    return widths,fmin,fmax
603   
604def getFWHM(TTh,Inst):
605    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
606    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
607    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
608    s = sig(TTh/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
609    g = gam(TTh/2.,Inst['X'][1],Inst['Y'][1])*100.
610    return gamFW(g,s)   
611               
612def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
613    DX = xdata[1]-xdata[0]
614    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
615    x = np.linspace(pos-fmin,pos+fmin,256)
616    dx = x[1]-x[0]
617    Norm = norm.pdf(x,loc=pos,scale=widths[0])
618    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
619    arg = [pos,shl/57.2958,dx,]
620    FCJ = fcjde.pdf(x,*arg,loc=pos)
621    if len(np.nonzero(FCJ)[0])>5:
622        z = np.column_stack([Norm,Cauchy,FCJ]).T
623        Z = fft(z)
624        Df = ifft(Z.prod(axis=0)).real
625    else:
626        z = np.column_stack([Norm,Cauchy]).T
627        Z = fft(z)
628        Df = fftshift(ifft(Z.prod(axis=0))).real
629    Df /= np.sum(Df)
630    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
631    return intens*Df(xdata)*DX/dx
632
633def getBackground(pfx,parmDict,bakType,xdata):
634    yb = np.zeros_like(xdata)
635    nBak = 0
636    while True:
637        key = pfx+'Back:'+str(nBak)
638        if key in parmDict:
639            nBak += 1
640        else:
641            break
642    if bakType in ['chebyschev','cosine']:
643        for iBak in range(nBak):   
644            key = pfx+'Back:'+str(iBak)
645            if bakType == 'chebyschev':
646                yb += parmDict[key]*(xdata-xdata[0])**iBak
647            elif bakType == 'cosine':
648                yb += parmDict[key]*npcosd(xdata*iBak)
649    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
650        if nBak == 1:
651            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
652        elif nBak == 2:
653            dX = xdata[-1]-xdata[0]
654            T2 = (xdata-xdata[0])/dX
655            T1 = 1.0-T2
656            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
657        else:
658            if bakType == 'lin interpolate':
659                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
660            elif bakType == 'inv interpolate':
661                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
662            elif bakType == 'log interpolate':
663                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
664            bakPos[0] = xdata[0]
665            bakPos[-1] = xdata[-1]
666            bakVals = np.zeros(nBak)
667            for i in range(nBak):
668                bakVals[i] = parmDict[pfx+'Back:'+str(i)]
669            bakInt = si.interp1d(bakPos,bakVals,'linear')
670            yb = bakInt(xdata)
671    if 'difC' in parmDict:
672        ff = 1.
673    else:       
674        try:
675            wave = parmDict[pfx+'Lam']
676        except KeyError:
677            wave = parmDict[pfx+'Lam1']
678        q = 4.0*np.pi*npsind(xdata/2.0)/wave
679        SQ = (q/(4.*np.pi))**2
680        FF = G2elem.GetFormFactorCoeff('Si')[0]
681        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
682    iD = 0       
683    while True:
684        try:
685            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
686            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
687            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
688            yb += ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
689            iD += 1       
690        except KeyError:
691            break
692    iD = 0
693    while True:
694        try:
695            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
696            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
697            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
698            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
699            shl = 0.002
700            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
701            iBeg = np.searchsorted(xdata,pkP-fmin)
702            iFin = np.searchsorted(xdata,pkP+fmax)
703            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
704            iD += 1       
705        except KeyError:
706            break       
707    return yb
708   
709def getBackgroundDerv(pfx,parmDict,bakType,xdata):
710    nBak = 0
711    while True:
712        key = pfx+'Back:'+str(nBak)
713        if key in parmDict:
714            nBak += 1
715        else:
716            break
717    dydb = np.zeros(shape=(nBak,len(xdata)))
718    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
719    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
720    dx = xdata[1]-xdata[0]
721
722    if bakType in ['chebyschev','cosine']:
723        for iBak in range(nBak):   
724            if bakType == 'chebyschev':
725                dydb[iBak] = (xdata-xdata[0])**iBak
726            elif bakType == 'cosine':
727                dydb[iBak] = npcosd(xdata*iBak)
728    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
729        if nBak == 1:
730            dydb[0] = np.ones_like(xdata)
731        elif nBak == 2:
732            dX = xdata[-1]-xdata[0]
733            T2 = (xdata-xdata[0])/dX
734            T1 = 1.0-T2
735            dydb = [T1,T2]
736        else:
737            if bakType == 'lin interpolate':
738                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
739            elif bakType == 'inv interpolate':
740                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
741            elif bakType == 'log interpolate':
742                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
743            bakPos[0] = xdata[0]
744            bakPos[-1] = xdata[-1]
745            dx = bakPos[1]-bakPos[0]
746            for i,pos in enumerate(bakPos):
747                if i == 0:
748                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
749                elif i == len(bakPos)-1:
750                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
751                else:
752                    dydb[i] = np.where(xdata>bakPos[i],
753                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
754                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
755    if 'difC' in parmDict:
756        ff = 1.
757    else:
758        try:
759            wave = parmDict[pfx+'Lam']
760        except KeyError:
761            wave = parmDict[pfx+'Lam1']
762        q = 4.0*np.pi*npsind(xdata/2.0)/wave
763        SQ = (q/(4*np.pi))**2
764        FF = G2elem.GetFormFactorCoeff('Si')[0]
765        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
766    iD = 0       
767    while True:
768        try:
769            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
770            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
771            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
772            sqr = np.sin(q*dbR)/(q*dbR)
773            cqr = np.cos(q*dbR)
774            temp = np.exp(-dbU*q**2)
775            dyddb[3*iD] = ff*sqr*temp/(np.pi*dx)
776            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*dx)
777            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*dx)
778            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
779        except KeyError:
780            break
781    iD = 0
782    while True:
783        try:
784            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
785            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
786            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
787            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
788            shl = 0.002
789            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
790            iBeg = np.searchsorted(xdata,pkP-fmin)
791            iFin = np.searchsorted(xdata,pkP+fmax)
792            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
793            dydpk[4*iD][iBeg:iFin] += 100.*dx*pkI*dFdp
794            dydpk[4*iD+1][iBeg:iFin] += 100.*dx*Df
795            dydpk[4*iD+2][iBeg:iFin] += 100.*dx*pkI*dFds
796            dydpk[4*iD+3][iBeg:iFin] += 100.*dx*pkI*dFdg
797            iD += 1       
798        except KeyError:
799            break
800    return dydb,dyddb,dydpk
801
802#use old fortran routine
803def getFCJVoigt3(pos,sig,gam,shl,xdata):
804   
805    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
806#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
807    Df /= np.sum(Df)
808    return Df
809
810def getdFCJVoigt3(pos,sig,gam,shl,xdata):
811   
812    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
813#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
814    sumDf = np.sum(Df)
815    return Df,dFdp,dFds,dFdg,dFdsh
816
817def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
818    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
819    Df /= np.sum(Df)
820    return Df 
821   
822def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
823    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
824    sumDf = np.sum(Df)
825    return Df,dFdp,dFda,dFdb,dFds,dFdg   
826
827def ellipseSize(H,Sij,GB):
828    HX = np.inner(H.T,GB)
829    lenHX = np.sqrt(np.sum(HX**2))
830    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
831    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
832    lenR = np.sqrt(np.sum(R**2))
833    return lenR
834
835def ellipseSizeDerv(H,Sij,GB):
836    lenR = ellipseSize(H,Sij,GB)
837    delt = 0.001
838    dRdS = np.zeros(6)
839    for i in range(6):
840        dSij = Sij[:]
841        dSij[i] += delt
842        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
843    return lenR,dRdS
844
845def getHKLpeak(dmin,SGData,A):
846    HKL = G2lat.GenHLaue(dmin,SGData,A)       
847    HKLs = []
848    for h,k,l,d in HKL:
849        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
850        if not ext:
851            HKLs.append([h,k,l,d,-1])
852    return HKLs
853
854def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
855   
856    yb = getBackground('',parmDict,bakType,xdata)
857    yc = np.zeros_like(yb)
858    if 'C' in dataType:
859        dx = xdata[1]-xdata[0]
860        U = parmDict['U']
861        V = parmDict['V']
862        W = parmDict['W']
863        X = parmDict['X']
864        Y = parmDict['Y']
865        shl = max(parmDict['SH/L'],0.002)
866        Ka2 = False
867        if 'Lam1' in parmDict.keys():
868            Ka2 = True
869            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
870            kRatio = parmDict['I(L2)/I(L1)']
871        iPeak = 0
872        while True:
873            try:
874                pos = parmDict['pos'+str(iPeak)]
875                theta = (pos-parmDict['Zero'])/2.0
876                intens = parmDict['int'+str(iPeak)]
877                sigName = 'sig'+str(iPeak)
878                if sigName in varyList:
879                    sig = parmDict[sigName]
880                else:
881                    sig = U*tand(theta)**2+V*tand(theta)+W
882                sig = max(sig,0.001)          #avoid neg sigma
883                gamName = 'gam'+str(iPeak)
884                if gamName in varyList:
885                    gam = parmDict[gamName]
886                else:
887                    gam = X/cosd(theta)+Y*tand(theta)
888                gam = max(gam,0.001)             #avoid neg gamma
889                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
890                iBeg = np.searchsorted(xdata,pos-fmin)
891                lenX = len(xdata)
892                if not iBeg:
893                    iFin = np.searchsorted(xdata,pos+fmin)
894                elif iBeg == lenX:
895                    iFin = iBeg
896                else:
897                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
898                if not iBeg+iFin:       #peak below low limit
899                    iPeak += 1
900                    continue
901                elif not iBeg-iFin:     #peak above high limit
902                    return yb+yc
903                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
904                if Ka2:
905                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
906                    kdelt = int((pos2-pos)/dx)               
907                    iBeg = min(lenX,iBeg+kdelt)
908                    iFin = min(lenX,iFin+kdelt)
909                    if iBeg-iFin:
910                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
911                iPeak += 1
912            except KeyError:        #no more peaks to process
913                return yb+yc
914    else:
915        difC = parmDict['difC']
916        alp0 = parmDict['alpha']
917        bet0 = parmDict['beta-0']
918        bet1 = parmDict['beta-1']
919        sig0 = parmDict['var-inst']
920        X = parmDict['X']
921        Y = parmDict['Y']
922        iPeak = 0
923        while True:
924            try:
925                pos = parmDict['pos'+str(iPeak)]               
926                tof = pos-parmDict['Zero']
927                dsp = tof/difC
928                intens = parmDict['int'+str(iPeak)]
929                alpName = 'alp'+str(iPeak)
930                if alpName in varyList:
931                    alp = parmDict[alpName]
932                else:
933                    alp = alp0/dsp
934                betName = 'bet'+str(iPeak)
935                if betName in varyList:
936                    bet = parmDict[betName]
937                else:
938                    bet = bet0+bet1/dsp**4
939                sigName = 'sig'+str(iPeak)
940                if sigName in varyList:
941                    sig = parmDict[sigName]
942                else:
943                    sig = sig0*dsp**2
944                gamName = 'gam'+str(iPeak)
945                if gamName in varyList:
946                    gam = parmDict[gamName]
947                else:
948                    gam = X*dsp**2+Y*dsp
949                gam = max(gam,0.001)             #avoid neg gamma
950                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
951                iBeg = np.searchsorted(xdata,pos-fmin)
952                iFin = np.searchsorted(xdata,pos+fmax)
953                lenX = len(xdata)
954                if not iBeg:
955                    iFin = np.searchsorted(xdata,pos+fmax)
956                elif iBeg == lenX:
957                    iFin = iBeg
958                else:
959                    iFin = np.searchsorted(xdata,pos+fmax)
960                if not iBeg+iFin:       #peak below low limit
961                    iPeak += 1
962                    continue
963                elif not iBeg-iFin:     #peak above high limit
964                    return yb+yc
965                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
966                iPeak += 1
967            except KeyError:        #no more peaks to process
968                return yb+yc
969           
970def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
971# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
972    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
973    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
974    if 'Back:0' in varyList:            #background derivs are in front if present
975        dMdv[0:len(dMdb)] = dMdb
976    names = ['DebyeA','DebyeR','DebyeU']
977    for name in varyList:
978        if 'Debye' in name:
979            parm,id = name.split(':')
980            ip = names.index(parm)
981            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
982    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
983    for name in varyList:
984        if 'BkPk' in name:
985            parm,id = name.split(':')
986            ip = names.index(parm)
987            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
988    if 'C' in dataType:
989        dx = xdata[1]-xdata[0]
990        U = parmDict['U']
991        V = parmDict['V']
992        W = parmDict['W']
993        X = parmDict['X']
994        Y = parmDict['Y']
995        shl = max(parmDict['SH/L'],0.002)
996        Ka2 = False
997        if 'Lam1' in parmDict.keys():
998            Ka2 = True
999            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1000            kRatio = parmDict['I(L2)/I(L1)']
1001        iPeak = 0
1002        while True:
1003            try:
1004                pos = parmDict['pos'+str(iPeak)]
1005                theta = (pos-parmDict['Zero'])/2.0
1006                intens = parmDict['int'+str(iPeak)]
1007                sigName = 'sig'+str(iPeak)
1008                tanth = tand(theta)
1009                costh = cosd(theta)
1010                if sigName in varyList:
1011                    sig = parmDict[sigName]
1012                else:
1013                    sig = U*tanth**2+V*tanth+W
1014                    dsdU = tanth**2
1015                    dsdV = tanth
1016                    dsdW = 1.0
1017                sig = max(sig,0.001)          #avoid neg sigma
1018                gamName = 'gam'+str(iPeak)
1019                if gamName in varyList:
1020                    gam = parmDict[gamName]
1021                else:
1022                    gam = X/costh+Y*tanth
1023                    dgdX = 1.0/costh
1024                    dgdY = tanth
1025                gam = max(gam,0.001)             #avoid neg gamma
1026                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1027                iBeg = np.searchsorted(xdata,pos-fmin)
1028                lenX = len(xdata)
1029                if not iBeg:
1030                    iFin = np.searchsorted(xdata,pos+fmin)
1031                elif iBeg == lenX:
1032                    iFin = iBeg
1033                else:
1034                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
1035                if not iBeg+iFin:       #peak below low limit
1036                    iPeak += 1
1037                    continue
1038                elif not iBeg-iFin:     #peak above high limit
1039                    break
1040                dMdpk = np.zeros(shape=(6,len(xdata)))
1041                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1042                for i in range(1,5):
1043                    dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i]
1044                dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0]
1045                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1046                if Ka2:
1047                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1048                    kdelt = int((pos2-pos)/dx)               
1049                    iBeg = min(lenX,iBeg+kdelt)
1050                    iFin = min(lenX,iFin+kdelt)
1051                    if iBeg-iFin:
1052                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1053                        for i in range(1,5):
1054                            dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i]
1055                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
1056                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
1057                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1058                for parmName in ['pos','int','sig','gam']:
1059                    try:
1060                        idx = varyList.index(parmName+str(iPeak))
1061                        dMdv[idx] = dervDict[parmName]
1062                    except ValueError:
1063                        pass
1064                if 'U' in varyList:
1065                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1066                if 'V' in varyList:
1067                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1068                if 'W' in varyList:
1069                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1070                if 'X' in varyList:
1071                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1072                if 'Y' in varyList:
1073                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1074                if 'SH/L' in varyList:
1075                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1076                if 'I(L2)/I(L1)' in varyList:
1077                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1078                iPeak += 1
1079            except KeyError:        #no more peaks to process
1080                break
1081    else:
1082        difC = parmDict['difC']
1083        alp0 = parmDict['alpha']
1084        bet0 = parmDict['beta-0']
1085        bet1 = parmDict['beta-1']
1086        sig0 = parmDict['var-inst']
1087        X = parmDict['X']
1088        Y = parmDict['Y']
1089        iPeak = 0
1090        while True:
1091            try:
1092                pos = parmDict['pos'+str(iPeak)]               
1093                tof = pos-parmDict['Zero']
1094                dsp = tof/difC
1095                intens = parmDict['int'+str(iPeak)]
1096                alpName = 'alp'+str(iPeak)
1097                if alpName in varyList:
1098                    alp = parmDict[alpName]
1099                else:
1100                    alp = alp0/dsp
1101                    dada0 = 1./dsp
1102                betName = 'bet'+str(iPeak)
1103                if betName in varyList:
1104                    bet = parmDict[betName]
1105                else:
1106                    bet = bet0+bet1/dsp**4
1107                    dbdb0 = 1.
1108                    dbdb1 = 1/dsp**4
1109                sigName = 'sig'+str(iPeak)
1110                if sigName in varyList:
1111                    sig = parmDict[sigName]
1112                else:
1113                    sig = sig0*dsp**2
1114                    dsds0 = dsp**2
1115                gamName = 'gam'+str(iPeak)
1116                if gamName in varyList:
1117                    gam = parmDict[gamName]
1118                else:
1119                    gam = X*dsp**2+Y*dsp
1120                    dsdX = dsp**2
1121                    dsdY = dsp
1122                gam = max(gam,0.001)             #avoid neg gamma
1123                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1124                iBeg = np.searchsorted(xdata,pos-fmin)
1125                lenX = len(xdata)
1126                if not iBeg:
1127                    iFin = np.searchsorted(xdata,pos+fmax)
1128                elif iBeg == lenX:
1129                    iFin = iBeg
1130                else:
1131                    iFin = np.searchsorted(xdata,pos+fmax)
1132                if not iBeg+iFin:       #peak below low limit
1133                    iPeak += 1
1134                    continue
1135                elif not iBeg-iFin:     #peak above high limit
1136                    break
1137                dMdpk = np.zeros(shape=(7,len(xdata)))
1138                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1139                for i in range(1,6):
1140                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
1141                dMdpk[0][iBeg:iFin] += dMdipk[0]
1142                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1143                for parmName in ['pos','int','alp','bet','sig','gam']:
1144                    try:
1145                        idx = varyList.index(parmName+str(iPeak))
1146                        dMdv[idx] = dervDict[parmName]
1147                    except ValueError:
1148                        pass
1149                if 'alpha' in varyList:
1150                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1151                if 'beta-0' in varyList:
1152                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1153                if 'beta-1' in varyList:
1154                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1155                if 'var-inst' in varyList:
1156                    dMdv[varyList.index('var-inst')] += dsds0*dervDict['sig']
1157                if 'X' in varyList:
1158                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1159                if 'Y' in varyList:
1160                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1161                iPeak += 1
1162            except KeyError:        #no more peaks to process
1163                break
1164    return dMdv
1165       
1166def Dict2Values(parmdict, varylist):
1167    '''Use before call to leastsq to setup list of values for the parameters
1168    in parmdict, as selected by key in varylist'''
1169    return [parmdict[key] for key in varylist] 
1170   
1171def Values2Dict(parmdict, varylist, values):
1172    ''' Use after call to leastsq to update the parameter dictionary with
1173    values corresponding to keys in varylist'''
1174    parmdict.update(zip(varylist,values))
1175   
1176def SetBackgroundParms(Background):
1177    if len(Background) == 1:            # fix up old backgrounds
1178        BackGround.append({'nDebye':0,'debyeTerms':[]})
1179    bakType,bakFlag = Background[0][:2]
1180    backVals = Background[0][3:]
1181    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1182    Debye = Background[1]           #also has background peaks stuff
1183    backDict = dict(zip(backNames,backVals))
1184    backVary = []
1185    if bakFlag:
1186        backVary = backNames
1187
1188    backDict['nDebye'] = Debye['nDebye']
1189    debyeDict = {}
1190    debyeList = []
1191    for i in range(Debye['nDebye']):
1192        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1193        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1194        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1195    debyeVary = []
1196    for item in debyeList:
1197        if item[1]:
1198            debyeVary.append(item[0])
1199    backDict.update(debyeDict)
1200    backVary += debyeVary
1201
1202    backDict['nPeaks'] = Debye['nPeaks']
1203    peaksDict = {}
1204    peaksList = []
1205    for i in range(Debye['nPeaks']):
1206        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1207        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1208        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1209    peaksVary = []
1210    for item in peaksList:
1211        if item[1]:
1212            peaksVary.append(item[0])
1213    backDict.update(peaksDict)
1214    backVary += peaksVary   
1215    return bakType,backDict,backVary
1216           
1217def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1218   
1219       
1220    def GetBackgroundParms(parmList,Background):
1221        iBak = 0
1222        while True:
1223            try:
1224                bakName = 'Back:'+str(iBak)
1225                Background[0][iBak+3] = parmList[bakName]
1226                iBak += 1
1227            except KeyError:
1228                break
1229        iDb = 0
1230        while True:
1231            names = ['DebyeA:','DebyeR:','DebyeU:']
1232            try:
1233                for i,name in enumerate(names):
1234                    val = parmList[name+str(iDb)]
1235                    Background[1]['debyeTerms'][iDb][2*i] = val
1236                iDb += 1
1237            except KeyError:
1238                break
1239        iDb = 0
1240        while True:
1241            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1242            try:
1243                for i,name in enumerate(names):
1244                    val = parmList[name+str(iDb)]
1245                    Background[1]['peaksList'][iDb][2*i] = val
1246                iDb += 1
1247            except KeyError:
1248                break
1249               
1250    def BackgroundPrint(Background,sigDict):
1251        if Background[0][1]:
1252            print 'Background coefficients for',Background[0][0],'function'
1253            ptfmt = "%12.5f"
1254            ptstr =  'values:'
1255            sigstr = 'esds  :'
1256            for i,back in enumerate(Background[0][3:]):
1257                ptstr += ptfmt % (back)
1258                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1259            print ptstr
1260            print sigstr
1261        else:
1262            print 'Background not refined'
1263        if Background[1]['nDebye']:
1264            parms = ['DebyeA','DebyeR','DebyeU']
1265            print 'Debye diffuse scattering coefficients'
1266            ptfmt = "%12.5f"
1267            names =   'names :'
1268            ptstr =  'values:'
1269            sigstr = 'esds  :'
1270            for item in sigDict:
1271                if 'Debye' in item:
1272                    names += '%12s'%(item)
1273                    sigstr += ptfmt%(sigDict[item])
1274                    parm,id = item.split(':')
1275                    ip = parms.index(parm)
1276                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1277            print names
1278            print ptstr
1279            print sigstr
1280        if Background[1]['nPeaks']:
1281            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1282            print 'Peaks in background coefficients'
1283            ptfmt = "%12.5f"
1284            names =   'names :'
1285            ptstr =  'values:'
1286            sigstr = 'esds  :'
1287            for item in sigDict:
1288                if 'BkPk' in item:
1289                    names += '%12s'%(item)
1290                    sigstr += ptfmt%(sigDict[item])
1291                    parm,id = item.split(':')
1292                    ip = parms.index(parm)
1293                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1294            print names
1295            print ptstr
1296            print sigstr
1297                           
1298    def SetInstParms(Inst):
1299        dataType = Inst['Type'][0]
1300        insVary = []
1301        insNames = []
1302        insVals = []
1303        for parm in Inst:
1304            insNames.append(parm)
1305            insVals.append(Inst[parm][1])
1306            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1307                'beta-0','beta-1','var-inst',] and Inst[parm][2]:
1308                    insVary.append(parm)
1309        instDict = dict(zip(insNames,insVals))
1310        instDict['X'] = max(instDict['X'],0.01)
1311        instDict['Y'] = max(instDict['Y'],0.01)
1312        if 'SH/L' in instDict:
1313            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1314        return dataType,instDict,insVary
1315       
1316    def GetInstParms(parmDict,Inst,varyList,Peaks):
1317        for name in Inst:
1318            Inst[name][1] = parmDict[name]
1319        iPeak = 0
1320        while True:
1321            try:
1322                sigName = 'sig'+str(iPeak)
1323                pos = parmDict['pos'+str(iPeak)]
1324                if sigName not in varyList:
1325                    if 'C' in Inst['Type'][0]:
1326                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1327                    else:
1328                        dsp = pos/Inst['difC'][1]
1329                        parmDict[sigName] = parmDict['var-inst']*dsp**2
1330                gamName = 'gam'+str(iPeak)
1331                if gamName not in varyList:
1332                    if 'C' in Inst['Type'][0]:
1333                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1334                    else:
1335                        dsp = pos/Inst['difC'][1]
1336                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1337                iPeak += 1
1338            except KeyError:
1339                break
1340       
1341    def InstPrint(Inst,sigDict):
1342        print 'Instrument Parameters:'
1343        ptfmt = "%12.6f"
1344        ptlbls = 'names :'
1345        ptstr =  'values:'
1346        sigstr = 'esds  :'
1347        for parm in Inst:
1348            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1349                'beta-0','beta-1','var-inst',]:
1350                ptlbls += "%s" % (parm.center(12))
1351                ptstr += ptfmt % (Inst[parm][1])
1352                if parm in sigDict:
1353                    sigstr += ptfmt % (sigDict[parm])
1354                else:
1355                    sigstr += 12*' '
1356        print ptlbls
1357        print ptstr
1358        print sigstr
1359
1360    def SetPeaksParms(dataType,Peaks):
1361        peakNames = []
1362        peakVary = []
1363        peakVals = []
1364        if 'C' in dataType:
1365            names = ['pos','int','sig','gam']
1366        else:
1367            names = ['pos','int','alpha','beta','sig','gam']
1368        for i,peak in enumerate(Peaks):
1369            for j,name in enumerate(names):
1370                peakVals.append(peak[2*j])
1371                parName = name+str(i)
1372                peakNames.append(parName)
1373                if peak[2*j+1]:
1374                    peakVary.append(parName)
1375        return dict(zip(peakNames,peakVals)),peakVary
1376               
1377    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1378        if 'C' in Inst['Type'][0]:
1379            names = ['pos','int','sig','gam']
1380        else:
1381            names = ['pos','int','alpha','beta','sig','gam']
1382        for i,peak in enumerate(Peaks):
1383            pos = parmDict['pos'+str(i)]
1384            if 'difC' in Inst:
1385                dsp = pos/Inst['difC'][1]
1386            for j in range(len(names)):
1387                parName = names[j]+str(i)
1388                if parName in varyList:
1389                    peak[2*j] = parmDict[parName]
1390                elif 'alpha' in parName:
1391                    peak[2*j] = parmDict['alpha']/dsp
1392                elif 'beta' in parName:
1393                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
1394                elif 'sig' in parName:
1395                    if 'C' in Inst['Type'][0]:
1396                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1397                    else:
1398                        peak[2*j] = parmDict['var-inst']*dsp**2
1399                elif 'gam' in parName:
1400                    if 'C' in Inst['Type'][0]:
1401                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1402                    else:
1403                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1404                       
1405    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1406        print 'Peak coefficients:'
1407        if 'C' in dataType:
1408            names = ['pos','int','sig','gam']
1409        else:
1410            names = ['pos','int','alpha','beta','sig','gam']           
1411        head = 13*' '
1412        for name in names:
1413            head += name.center(10)+'esd'.center(10)
1414        print head
1415        if 'C' in dataType:
1416            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1417        else:
1418            ptfmt = {'pos':"%10.2f",'int':"%10.1f",'alpha':"%10.3f",'beta':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1419        for i,peak in enumerate(Peaks):
1420            ptstr =  ':'
1421            for j in range(len(names)):
1422                name = names[j]
1423                parName = name+str(i)
1424                ptstr += ptfmt[name] % (parmDict[parName])
1425                if parName in varyList:
1426#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1427                    ptstr += ptfmt[name] % (sigDict[parName])
1428                else:
1429#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1430                    ptstr += 10*' '
1431            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1432               
1433    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1434        parmdict.update(zip(varylist,values))
1435        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1436           
1437    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1438        parmdict.update(zip(varylist,values))
1439        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1440        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1441        if dlg:
1442            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1443            if not GoOn:
1444                return -M           #abort!!
1445        return M
1446       
1447    if controls:
1448        Ftol = controls['min dM/M']
1449        derivType = controls['deriv type']
1450    else:
1451        Ftol = 0.0001
1452        derivType = 'analytic'
1453    if oneCycle:
1454        Ftol = 1.0
1455    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1456    xBeg = np.searchsorted(x,Limits[0])
1457    xFin = np.searchsorted(x,Limits[1])
1458    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1459    dataType,insDict,insVary = SetInstParms(Inst)
1460    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1461    parmDict = {}
1462    parmDict.update(bakDict)
1463    parmDict.update(insDict)
1464    parmDict.update(peakDict)
1465    varyList = bakVary+insVary+peakVary
1466    while True:
1467        begin = time.time()
1468        values =  np.array(Dict2Values(parmDict, varyList))
1469        if FitPgm == 'LSQ':
1470            try:
1471                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1472                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1473                ncyc = int(result[2]['nfev']/2)
1474            finally:
1475                dlg.Destroy()
1476            runtime = time.time()-begin   
1477            chisq = np.sum(result[2]['fvec']**2)
1478            Values2Dict(parmDict, varyList, result[0])
1479            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1480            GOF = chisq/(xFin-xBeg-len(varyList))
1481            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1482            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1483            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1484            try:
1485                sig = np.sqrt(np.diag(result[1])*GOF)
1486                if np.any(np.isnan(sig)):
1487                    print '*** Least squares aborted - some invalid esds possible ***'
1488                break                   #refinement succeeded - finish up!
1489            except ValueError:          #result[1] is None on singular matrix
1490                print '**** Refinement failed - singular matrix ****'
1491                Ipvt = result[2]['ipvt']
1492                for i,ipvt in enumerate(Ipvt):
1493                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1494                        print 'Removing parameter: ',varyList[ipvt-1]
1495                        del(varyList[ipvt-1])
1496                        break
1497        elif FitPgm == 'BFGS':
1498            print 'Other program here'
1499            return
1500       
1501    sigDict = dict(zip(varyList,sig))
1502    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1503    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1504    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1505    GetBackgroundParms(parmDict,Background)
1506    BackgroundPrint(Background,sigDict)
1507    GetInstParms(parmDict,Inst,varyList,Peaks)
1508    InstPrint(Inst,sigDict)
1509    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1510    PeaksPrint(dataType,parmDict,sigDict,varyList)
1511
1512def calcIncident(Iparm,xdata):
1513
1514    def IfunAdv(Iparm,xdata):
1515        Itype = Iparm['Itype']
1516        Icoef = Iparm['Icoeff']
1517        DYI = np.ones((12,xdata.shape[0]))
1518        YI = np.ones_like(xdata)*Icoef[0]
1519       
1520        x = xdata/1000.                 #expressions are in ms
1521        if Itype == 'Exponential':
1522            for i in range(1,10,2):
1523                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1524                YI += Icoef[i]*Eterm
1525                DYI[i] *= Eterm
1526                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1527        elif 'Maxwell'in Itype:
1528            Eterm = np.exp(-Icoef[2]/x**2)
1529            DYI[1] = Eterm/x**5
1530            DYI[2] = -Icoef[1]*DYI[1]/x**2
1531            YI += (Icoef[1]*Eterm/x**5)
1532            if 'Exponential' in Itype:
1533                for i in range(3,12,2):
1534                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1535                    YI += Icoef[i]*Eterm
1536                    DYI[i] *= Eterm
1537                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1538            else:   #Chebyschev
1539                T = (2./x)-1.
1540                Ccof = np.ones((12,xdata.shape[0]))
1541                Ccof[1] = T
1542                for i in range(2,12):
1543                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1544                for i in range(1,10):
1545                    YI += Ccof[i]*Icoef[i+2]
1546                    DYI[i+2] =Ccof[i]
1547        return YI,DYI
1548       
1549    Iesd = np.array(Iparm['Iesd'])
1550    Icovar = Iparm['Icovar']
1551    YI,DYI = IfunAdv(Iparm,xdata)
1552    YI = np.where(YI>0,YI,1.)
1553    WYI = np.zeros_like(xdata)
1554    vcov = np.zeros((12,12))
1555    k = 0
1556    for i in range(12):
1557        for j in range(i,12):
1558            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1559            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1560            k += 1
1561    M = np.inner(vcov,DYI.T)
1562    WYI = np.sum(M*DYI,axis=0)
1563    WYI = np.where(WYI>0.,WYI,0.)
1564    return YI,WYI
1565   
1566#testing data
1567NeedTestData = True
1568def TestData():
1569#    global NeedTestData
1570    NeedTestData = False
1571    global bakType
1572    bakType = 'chebyschev'
1573    global xdata
1574    xdata = np.linspace(4.0,40.0,36000)
1575    global parmDict0
1576    parmDict0 = {
1577        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1578        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1579        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1580        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1581        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1582        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1583        }
1584    global parmDict1
1585    parmDict1 = {
1586        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1587        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1588        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1589        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1590        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1591        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1592        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1593        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1594        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1595        }
1596    global parmDict2
1597    parmDict2 = {
1598        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1599        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1600        'Back0':5.,'Back1':-0.02,'Back2':.004,
1601#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1602        }
1603    global varyList
1604    varyList = []
1605
1606def test0():
1607    if NeedTestData: TestData()
1608    msg = 'test '
1609    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1610    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1611    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1612    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1613    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1614    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1615   
1616def test1():
1617    if NeedTestData: TestData()
1618    time0 = time.time()
1619    for i in range(100):
1620        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1621    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1622   
1623def test2(name,delt):
1624    if NeedTestData: TestData()
1625    varyList = [name,]
1626    xdata = np.linspace(5.6,5.8,400)
1627    hplot = plotter.add('derivatives test for '+name).gca()
1628    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1629    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1630    parmDict2[name] += delt
1631    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1632    hplot.plot(xdata,(y1-y0)/delt,'r+')
1633   
1634def test3(name,delt):
1635    if NeedTestData: TestData()
1636    names = ['pos','sig','gam','shl']
1637    idx = names.index(name)
1638    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1639    xdata = np.linspace(5.6,5.8,800)
1640    dx = xdata[1]-xdata[0]
1641    hplot = plotter.add('derivatives test for '+name).gca()
1642    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1643    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1644    myDict[name] += delt
1645    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1646    hplot.plot(xdata,(y1-y0)/delt,'r+')
1647
1648if __name__ == '__main__':
1649    import GSASIItestplot as plot
1650    global plotter
1651    plotter = plot.PlotNotebook()
1652#    test0()
1653#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1654    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1655        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1656        test2(name,shft)
1657    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1658        test3(name,shft)
1659    print "OK"
1660    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.