source: branch/MPbranch/GSASIIpwd.py

Last change on this file was 507, checked in by toby, 10 years ago

fix function calc to allow for out of sequence finish of 'threads'; other minor fixes

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