source: trunk/GSASIIpwd.py @ 358

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

small fixes
File size: 44.5 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2011-04-20 13:09:53 -0500 (Wed, 20 Apr 2011) $
6# $Author: vondreele $
7# $Revision: 267 $
8# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 267 2011-04-20 18:09:53Z vondreele $
10########### SVN repository information ###################
11import sys
12import math
13import wx
14import time
15
16import numpy as np
17import scipy as sp
18import numpy.linalg as nl
19from numpy.fft import ifft, fft, fftshift
20import scipy.interpolate as si
21import scipy.stats as st
22import scipy.optimize as so
23
24import GSASIIpath
25import GSASIIplot as G2plt
26import GSASIIlattice as G2lat
27import GSASIIElem as G2elem
28import GSASIIgrid as G2gd
29import GSASIIIO as G2IO
30import pypowder as pyd
31
32# trig functions in degrees
33sind = lambda x: math.sin(x*math.pi/180.)
34asind = lambda x: 180.*math.asin(x)/math.pi
35tand = lambda x: math.tan(x*math.pi/180.)
36atand = lambda x: 180.*math.atan(x)/math.pi
37atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
38cosd = lambda x: math.cos(x*math.pi/180.)
39acosd = lambda x: 180.*math.acos(x)/math.pi
40rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
41#numpy versions
42npsind = lambda x: np.sin(x*np.pi/180.)
43npasind = lambda x: 180.*np.arcsin(x)/math.pi
44npcosd = lambda x: np.cos(x*math.pi/180.)
45npacosd = lambda x: 180.*np.arccos(x)/math.pi
46nptand = lambda x: np.tan(x*math.pi/180.)
47npatand = lambda x: 180.*np.arctan(x)/np.pi
48npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
49npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
50npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
51   
52#GSASII pdf calculation routines
53       
54def Transmission(Geometry,Abs,Diam):
55#Calculate sample transmission
56#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
57#   Abs: absorption coeff in cm-1
58#   Diam: sample thickness/diameter in mm
59    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
60        MuR = Abs*Diam/20.0
61        if MuR <= 3.0:
62            T0 = 16/(3.*math.pi)
63            T1 = -0.045780
64            T2 = -0.02489
65            T3 = 0.003045
66            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
67            if T < -20.:
68                return 2.06e-9
69            else:
70                return math.exp(T)
71        else:
72            T1 = 1.433902
73            T2 = 0.013869+0.337894
74            T3 = 1.933433+1.163198
75            T4 = 0.044365-0.04259
76            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
77            return T/100.
78    elif 'plate' in Geometry:
79        MuR = Abs*Diam/10.
80        return math.exp(-MuR)
81    elif 'Bragg' in Geometry:
82        return 0.0
83
84def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
85#Calculate sample absorption
86#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
87#   Abs: absorption coeff in cm-1
88#   Diam: sample thickness/diameter in mm
89#   Tth: 2-theta scattering angle - can be numpy array
90#   Phi: flat plate tilt angle - future
91#   Psi: flat plate tilt axis - future
92    Sth2 = npsind(Tth/2.0)**2
93    Cth2 = 1.-Sth2
94    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
95        MuR = Abs*Diam/20.0
96        if MuR < 3.0:
97            T0 = 16.0/(3*np.pi)
98            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
99                0.109561*np.sqrt(Sth2)-26.04556
100            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
101                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
102            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
103            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
104            return np.exp(Trns)
105        else:
106            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
107                10.02088*Sth2**3-3.36778*Sth2**4
108            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
109                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
110            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
111                0.13576*np.sqrt(Sth2)+1.163198
112            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
113            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
114            return Trns/100.
115    elif 'Bragg' in Geometry:
116        return 1.0
117    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
118        # and only defined for 2theta < 90
119        MuR = Abs*Diam/10.0
120        T1 = np.exp(-MuR)
121        T2 = np.exp(-MuR/npcosd(Tth))
122        Tb = MuR-MuR/npcosd(Tth)
123        return (T2-T1)/Tb
124    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
125        MuR = Abs*Diam/10.0
126        cth = npcosd(Tth/2.0)
127        return np.exp(-MuR/cth)/cth
128       
129def Polarization(Pola,Tth,Azm=0.0):
130#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
131#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
132#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
133#   Tth: 2-theta scattering angle - can be numpy array
134#       which (if either) of these is "right"?
135#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
136#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
137    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
138        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
139    dpdPola = npsind(Azm)**2*npsind(Tth)**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    if bakType == 'chebyschev':
610        iBak = 0
611        while True:
612            key = pfx+'Back:'+str(iBak)
613            try:
614                yb += parmDict[key]*(xdata-xdata[0])**iBak
615                iBak += 1
616            except KeyError:
617                break
618    return yb
619   
620def getBackgroundDerv(pfx,parmDict,bakType,xdata):
621    dydb = []
622    if bakType == 'chebyschev':
623        iBak = 0
624        while True:
625            if pfx+'Back:'+str(iBak) in parmDict:
626                dydb.append((xdata-xdata[0])**iBak)
627                iBak += 1
628            else:
629                break
630    return dydb
631
632#use old fortran routine
633def getFCJVoigt3(pos,sig,gam,shl,xdata):
634   
635    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
636    Df /= np.sum(Df)
637    return Df
638
639def getdFCJVoigt3(pos,sig,gam,shl,xdata):
640   
641    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
642    sumDf = np.sum(Df)
643    return Df,dFdp,dFds,dFdg,dFdsh
644   
645
646def getPeakProfile(parmDict,xdata,varyList,bakType):
647   
648    yb = getBackground('',parmDict,bakType,xdata)
649    yc = np.zeros_like(yb)
650    dx = xdata[1]-xdata[0]
651    U = parmDict['U']
652    V = parmDict['V']
653    W = parmDict['W']
654    X = parmDict['X']
655    Y = parmDict['Y']
656    shl = max(parmDict['SH/L'],0.002)
657    Ka2 = False
658    if 'Lam1' in parmDict.keys():
659        Ka2 = True
660        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
661        kRatio = parmDict['I(L2)/I(L1)']
662    iPeak = 0
663    while True:
664        try:
665            pos = parmDict['pos'+str(iPeak)]
666            intens = parmDict['int'+str(iPeak)]
667            sigName = 'sig'+str(iPeak)
668            if sigName in varyList:
669                sig = parmDict[sigName]
670            else:
671                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
672            sig = max(sig,0.001)          #avoid neg sigma
673            gamName = 'gam'+str(iPeak)
674            if gamName in varyList:
675                gam = parmDict[gamName]
676            else:
677                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
678            gam = max(gam,0.001)             #avoid neg gamma
679            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
680            iBeg = np.searchsorted(xdata,pos-fmin)
681            lenX = len(xdata)
682            if not iBeg:
683                iFin = np.searchsorted(xdata,pos+fmin)
684            elif iBeg == lenX:
685                iFin = iBeg
686            else:
687                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
688            if not iBeg+iFin:       #peak below low limit
689                iPeak += 1
690                continue
691            elif not iBeg-iFin:     #peak above high limit
692                return yb+yc
693            yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
694            if Ka2:
695                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
696                kdelt = int((pos2-pos)/dx)               
697                iBeg = min(lenX,iBeg+kdelt)
698                iFin = min(lenX,iFin+kdelt)
699                if iBeg-iFin:
700                    yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
701            iPeak += 1
702        except KeyError:        #no more peaks to process
703            return yb+yc
704           
705def getPeakProfileDerv(parmDict,xdata,varyList,bakType):
706# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
707    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
708    if 'Back:0' in varyList:            #background derivs are in front if present
709        dMdb = getBackgroundDerv('',parmDict,bakType,xdata)
710        dMdv[0:len(dMdb)] = dMdb
711       
712    dx = xdata[1]-xdata[0]
713    U = parmDict['U']
714    V = parmDict['V']
715    W = parmDict['W']
716    X = parmDict['X']
717    Y = parmDict['Y']
718    shl = max(parmDict['SH/L'],0.002)
719    Ka2 = False
720    if 'Lam1' in parmDict.keys():
721        Ka2 = True
722        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
723        kRatio = parmDict['I(L2)/I(L1)']
724    iPeak = 0
725    while True:
726        try:
727            pos = parmDict['pos'+str(iPeak)]
728            intens = parmDict['int'+str(iPeak)]
729            sigName = 'sig'+str(iPeak)
730            tanth = tand(pos/2.0)
731            costh = cosd(pos/2.0)
732            if sigName in varyList:
733                sig = parmDict[sigName]
734            else:
735                sig = U*tanth**2+V*tanth+W
736                dsdU = tanth**2
737                dsdV = tanth
738                dsdW = 1.0
739            sig = max(sig,0.001)          #avoid neg sigma
740            gamName = 'gam'+str(iPeak)
741            if gamName in varyList:
742                gam = parmDict[gamName]
743            else:
744                gam = X/costh+Y*tanth
745                dgdX = 1.0/costh
746                dgdY = tanth
747            gam = max(gam,0.001)             #avoid neg gamma
748            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
749            iBeg = np.searchsorted(xdata,pos-fmin)
750            lenX = len(xdata)
751            if not iBeg:
752                iFin = np.searchsorted(xdata,pos+fmin)
753            elif iBeg == lenX:
754                iFin = iBeg
755            else:
756                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
757            if not iBeg+iFin:       #peak below low limit
758                iPeak += 1
759                continue
760            elif not iBeg-iFin:     #peak above high limit
761                break
762            dMdpk = np.zeros(shape=(6,len(xdata)))
763            dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
764            for i in range(1,5):
765                dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i]
766            dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0]
767            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
768            if Ka2:
769                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
770                kdelt = int((pos2-pos)/dx)               
771                iBeg = min(lenX,iBeg+kdelt)
772                iFin = min(lenX,iFin+kdelt)
773                if iBeg-iFin:
774                    dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
775                    for i in range(1,5):
776                        dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i]
777                    dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
778                    dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
779                    dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
780            for parmName in ['pos','int','sig','gam']:
781                try:
782                    idx = varyList.index(parmName+str(iPeak))
783                    dMdv[idx] = dervDict[parmName]
784                except ValueError:
785                    pass
786            if 'U' in varyList:
787                dMdv[varyList.index('U')] += dsdU*dervDict['sig']
788            if 'V' in varyList:
789                dMdv[varyList.index('V')] += dsdV*dervDict['sig']
790            if 'W' in varyList:
791                dMdv[varyList.index('W')] += dsdW*dervDict['sig']
792            if 'X' in varyList:
793                dMdv[varyList.index('X')] += dgdX*dervDict['gam']
794            if 'Y' in varyList:
795                dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
796            if 'SH/L' in varyList:
797                dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
798            if 'I(L2)/I(L1)' in varyList:
799                dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
800            iPeak += 1
801        except KeyError:        #no more peaks to process
802            break
803    return dMdv
804       
805def Dict2Values(parmdict, varylist):
806    '''Use before call to leastsq to setup list of values for the parameters
807    in parmdict, as selected by key in varylist'''
808    return [parmdict[key] for key in varylist] 
809   
810def Values2Dict(parmdict, varylist, values):
811    ''' Use after call to leastsq to update the parameter dictionary with
812    values corresponding to keys in varylist'''
813    parmdict.update(zip(varylist,values))
814   
815def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False,controls=None):
816   
817    def SetBackgroundParms(Background):
818        bakType,bakFlag = Background[:2]
819        backVals = Background[3:]
820        backNames = ['Back:'+str(i) for i in range(len(backVals))]
821        if bakFlag: #returns backNames as varyList = backNames
822            return bakType,dict(zip(backNames,backVals)),backNames
823        else:       #no background varied; varyList = []
824            return bakType,dict(zip(backNames,backVals)),[]
825       
826    def GetBackgroundParms(parmList,Background):
827        iBak = 0
828        while True:
829            try:
830                bakName = 'Back:'+str(iBak)
831                Background[iBak+3] = parmList[bakName]
832                iBak += 1
833            except KeyError:
834                break
835               
836    def BackgroundPrint(Background,sigDict):
837        if Background[1]:
838            print 'Background coefficients for',Background[0],'function'
839            ptfmt = "%12.5f"
840            ptstr =  'values:'
841            sigstr = 'esds  :'
842            for i,back in enumerate(Background[3:]):
843                ptstr += ptfmt % (back)
844                sigstr += ptfmt % (sigDict['Back:'+str(i)])
845            print ptstr
846            print sigstr
847        else:
848            print 'Background not refined'
849           
850    def SetInstParms(Inst):
851        insVals,insFlags,insNames = Inst[1:4]
852        dataType = insVals[0]
853        insVary = []
854        for i,flag in enumerate(insFlags):
855            if flag and insNames[i] in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
856                insVary.append(insNames[i])
857        instDict = dict(zip(insNames,insVals))
858        instDict['X'] = max(instDict['X'],0.01)
859        instDict['Y'] = max(instDict['Y'],0.01)
860        instDict['SH/L'] = max(instDict['SH/L'],0.002)
861        return dataType,instDict,insVary
862       
863    def GetInstParms(parmDict,Inst,varyList,Peaks):
864        instNames = Inst[3]
865        for i,name in enumerate(instNames):
866            Inst[1][i] = parmDict[name]
867        iPeak = 0
868        while True:
869            try:
870                sigName = 'sig'+str(iPeak)
871                pos = parmDict['pos'+str(iPeak)]
872                if sigName not in varyList:
873                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
874                gamName = 'gam'+str(iPeak)
875                if gamName not in varyList:
876                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
877                iPeak += 1
878            except KeyError:
879                break
880       
881    def InstPrint(Inst,sigDict):
882        print 'Instrument Parameters:'
883        ptfmt = "%12.6f"
884        ptlbls = 'names :'
885        ptstr =  'values:'
886        sigstr = 'esds  :'
887        instNames = Inst[3][1:]
888        for i,name in enumerate(instNames):
889            ptlbls += "%s" % (name.center(12))
890            ptstr += ptfmt % (Inst[1][i+1])
891            if name in sigDict:
892                sigstr += ptfmt % (sigDict[name])
893            else:
894                sigstr += 12*' '
895        print ptlbls
896        print ptstr
897        print sigstr
898
899    def SetPeaksParms(Peaks):
900        peakNames = []
901        peakVary = []
902        peakVals = []
903        names = ['pos','int','sig','gam']
904        for i,peak in enumerate(Peaks):
905            for j in range(4):
906                peakVals.append(peak[2*j])
907                parName = names[j]+str(i)
908                peakNames.append(parName)
909                if peak[2*j+1]:
910                    peakVary.append(parName)
911        return dict(zip(peakNames,peakVals)),peakVary
912               
913    def GetPeaksParms(parmDict,Peaks,varyList):
914        names = ['pos','int','sig','gam']
915        for i,peak in enumerate(Peaks):
916            for j in range(4):
917                pos = parmDict['pos'+str(i)]
918                parName = names[j]+str(i)
919                if parName in varyList:
920                    peak[2*j] = parmDict[parName]
921                elif 'sig' in parName:
922                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
923                elif 'gam' in parName:
924                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
925                       
926    def PeaksPrint(parmDict,sigDict,varyList):
927        print 'Peak coefficients:'
928        names = ['pos','int','sig','gam']
929        head = 15*' '
930        for name in names:
931            head += name.center(12)+'esd'.center(12)
932        print head
933        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
934        for i,peak in enumerate(Peaks):
935            ptstr =  ':'
936            for j in range(4):
937                name = names[j]
938                parName = name+str(i)
939                ptstr += ptfmt[name] % (parmDict[parName])
940                if parName in varyList:
941#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
942                    ptstr += ptfmt[name] % (sigDict[parName])
943                else:
944#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
945                    ptstr += 12*' '
946            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
947               
948    def devPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):
949        parmdict.update(zip(varylist,values))
950        return np.sqrt(weights)*getPeakProfileDerv(parmdict,xdata,varylist,bakType)
951           
952    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
953        parmdict.update(zip(varylist,values))
954        M = np.sqrt(weights)*(getPeakProfile(parmdict,xdata,varylist,bakType)-ydata)
955        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
956        if dlg:
957            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
958            if not GoOn:
959                return -M           #abort!!
960        return M
961       
962    if controls:
963        Ftol = controls['min dM/M']
964        derivType = controls['deriv type']
965    else:
966        Ftol = 0.0001
967        derivType = 'analytic'
968    if oneCycle:
969        Ftol = 1.0
970    x,y,w,yc,yb,yd = data               #these are numpy arrays!
971    xBeg = np.searchsorted(x,Limits[0])
972    xFin = np.searchsorted(x,Limits[1])
973    bakType,bakDict,bakVary = SetBackgroundParms(Background)
974    dataType,insDict,insVary = SetInstParms(Inst)
975    peakDict,peakVary = SetPeaksParms(Peaks)
976    parmDict = {}
977    parmDict.update(bakDict)
978    parmDict.update(insDict)
979    parmDict.update(peakDict)
980    varyList = bakVary+insVary+peakVary
981    while True:
982        begin = time.time()
983        values =  np.array(Dict2Values(parmDict, varyList))
984        if FitPgm == 'LSQ':
985            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
986            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
987            screenSize = wx.ClientDisplayRect()
988            Size = dlg.GetSize()
989            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
990            try:
991                if derivType == 'analytic':
992                    result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
993                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
994                    ncyc = int(result[2]['nfev']/2)
995                else:
996                    result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,
997                        args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
998                    ncyc = int(result[2]['nfev']/len(varyList))
999            finally:
1000                dlg.Destroy()
1001            runtime = time.time()-begin   
1002            chisq = np.sum(result[2]['fvec']**2)
1003            Values2Dict(parmDict, varyList, result[0])
1004            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1005            GOF = chisq/(xFin-xBeg-len(varyList))
1006            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1007            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1008            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1009            try:
1010                sig = np.sqrt(np.diag(result[1])*GOF)
1011                if np.any(np.isnan(sig)):
1012                    print '*** Least squares aborted - some invalid esds possible ***'
1013                break                   #refinement succeeded - finish up!
1014            except ValueError:          #result[1] is None on singular matrix
1015                print '**** Refinement failed - singular matrix ****'
1016                Ipvt = result[2]['ipvt']
1017                for i,ipvt in enumerate(Ipvt):
1018                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1019                        print 'Removing parameter: ',varyList[ipvt-1]
1020                        del(varyList[ipvt-1])
1021                        break
1022        elif FitPgm == 'BFGS':
1023            print 'Other program here'
1024            return
1025       
1026    sigDict = dict(zip(varyList,sig))
1027    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1028    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
1029    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1030    GetBackgroundParms(parmDict,Background)
1031    BackgroundPrint(Background,sigDict)
1032    GetInstParms(parmDict,Inst,varyList,Peaks)
1033    InstPrint(Inst,sigDict)
1034    GetPeaksParms(parmDict,Peaks,varyList)   
1035    PeaksPrint(parmDict,sigDict,varyList)
1036   
1037#testing data
1038NeedTestData = True
1039def TestData():
1040#    global NeedTestData
1041    NeedTestData = False
1042    global bakType
1043    bakType = 'chebyschev'
1044    global xdata
1045    xdata = np.linspace(4.0,40.0,36000)
1046    global parmDict0
1047    parmDict0 = {
1048        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1049        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1050        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1051        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1052        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1053        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1054        }
1055    global parmDict1
1056    parmDict1 = {
1057        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1058        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1059        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1060        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1061        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1062        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1063        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1064        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1065        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1066        }
1067    global parmDict2
1068    parmDict2 = {
1069        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1070        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1071        'Back0':5.,'Back1':-0.02,'Back2':.004,
1072#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1073        }
1074    global varyList
1075    varyList = []
1076
1077def test0():
1078    if NeedTestData: TestData()
1079    msg = 'test '
1080    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1081    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1082    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1083    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1084    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1085    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1086   
1087def test1():
1088    if NeedTestData: TestData()
1089    time0 = time.time()
1090    for i in range(100):
1091        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1092    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1093   
1094def test2(name,delt):
1095    if NeedTestData: TestData()
1096    varyList = [name,]
1097    xdata = np.linspace(5.6,5.8,400)
1098    hplot = plotter.add('derivatives test for '+name).gca()
1099    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1100    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1101    parmDict2[name] += delt
1102    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1103    hplot.plot(xdata,(y1-y0)/delt,'r+')
1104   
1105def test3(name,delt):
1106    if NeedTestData: TestData()
1107    names = ['pos','sig','gam','shl']
1108    idx = names.index(name)
1109    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1110    xdata = np.linspace(5.6,5.8,800)
1111    dx = xdata[1]-xdata[0]
1112    hplot = plotter.add('derivatives test for '+name).gca()
1113    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1114    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1115    myDict[name] += delt
1116    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1117    hplot.plot(xdata,(y1-y0)/delt,'r+')
1118
1119if __name__ == '__main__':
1120    import GSASIItestplot as plot
1121    global plotter
1122    plotter = plot.PlotNotebook()
1123#    test0()
1124#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1125    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1126        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1127        test2(name,shft)
1128    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1129        test3(name,shft)
1130    print "OK"
1131    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.