source: trunk/GSASIIpwd.py @ 453

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

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