source: trunk/GSASIIpwd.py @ 432

Last change on this file since 432 was 432, checked in by vondreele, 11 years ago

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