source: trunk/GSASIIpwd.py @ 484

Last change on this file since 484 was 482, checked in by vondreele, 13 years ago

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