source: trunk/GSASIIpwd.py @ 350

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

implement fortran psVoigt routine - faster by 2X over python fft version

File size: 36.9 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       
53def factorize(num):
54    ''' Provide prime number factors for integer num
55    Returns dictionary of prime factors (keys) & power for each (data)
56    '''
57    factors = {}
58    orig = num
59
60    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
61    i, sqi = 2, 4
62    while sqi <= num:
63        while not num%i:
64            num /= i
65            factors[i] = factors.get(i, 0) + 1
66
67        sqi += 2*i + 1
68        i += 1
69
70    if num != 1 and num != orig:
71        factors[num] = factors.get(num, 0) + 1
72
73    if factors:
74        return factors
75    else:
76        return {num:1}          #a prime number!
77           
78def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
79    ''' Provide list of optimal data sizes for FFT calculations
80    Input:
81        nmin: minimum data size >= 1
82        nmax: maximum data size > nmin
83        thresh: maximum prime factor allowed
84    Returns:
85        list of data sizes where the maximum prime factor is < thresh
86    ''' 
87    plist = []
88    nmin = max(1,nmin)
89    nmax = max(nmin+1,nmax)
90    for p in range(nmin,nmax):
91        if max(factorize(p).keys()) < thresh:
92            plist.append(p)
93    return plist
94
95def Transmission(Geometry,Abs,Diam):
96#Calculate sample transmission
97#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
98#   Abs: absorption coeff in cm-1
99#   Diam: sample thickness/diameter in mm
100    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
101        MuR = Abs*Diam/20.0
102        if MuR <= 3.0:
103            T0 = 16/(3.*math.pi)
104            T1 = -0.045780
105            T2 = -0.02489
106            T3 = 0.003045
107            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
108            if T < -20.:
109                return 2.06e-9
110            else:
111                return math.exp(T)
112        else:
113            T1 = 1.433902
114            T2 = 0.013869+0.337894
115            T3 = 1.933433+1.163198
116            T4 = 0.044365-0.04259
117            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
118            return T/100.
119    elif 'plate' in Geometry:
120        MuR = Abs*Diam/10.
121        return math.exp(-MuR)
122    elif 'Bragg' in Geometry:
123        return 0.0
124
125def Absorb(Geometry,Abs,Diam,Tth,Phi=0,Psi=0):
126#Calculate sample absorption
127#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
128#   Abs: absorption coeff in cm-1
129#   Diam: sample thickness/diameter in mm
130#   Tth: 2-theta scattering angle - can be numpy array
131#   Phi: flat plate tilt angle - future
132#   Psi: flat plate tilt axis - future
133    Sth2 = npsind(Tth/2.0)**2
134    Cth2 = 1.-Sth2
135    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
136        MuR = Abs*Diam/20.0
137        if MuR < 3.0:
138            T0 = 16.0/(3*np.pi)
139            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
140                0.109561*np.sqrt(Sth2)-26.04556
141            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
142                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
143            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
144            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
145            return np.exp(Trns)
146        else:
147            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
148                10.02088*Sth2**3-3.36778*Sth2**4
149            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
150                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
151            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
152                0.13576*np.sqrt(Sth2)+1.163198
153            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
154            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
155            return Trns/100.
156    elif 'Bragg' in Geometry:
157        return 1.0
158    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
159        # and only defined for 2theta < 90
160        MuR = Abs*Diam/10.0
161        T1 = np.exp(-MuR)
162        T2 = np.exp(-MuR/npcosd(Tth))
163        Tb = MuR-MuR/npcosd(Tth)
164        return (T2-T1)/Tb
165    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
166        MuR = Abs*Diam/10.0
167        cth = npcosd(Tth/2.0)
168        return np.exp(-MuR/cth)/cth
169       
170def Polarization(Pola,Tth,Azm=0.0):
171#   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
172#   Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
173#   Azm: azimuthal angle e.g. 0.0 in plane of polarization
174#   Tth: 2-theta scattering angle - can be numpy array
175#       which (if either) of these is "right"?
176#    return (Pola*npcosd(Azm)**2+(1.-Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \
177#        Pola*npsind(Azm)**2+(1.-Pola)*npcosd(Azm)**2
178    return ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
179        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
180   
181def Oblique(ObCoeff,Tth):
182    if ObCoeff:
183        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
184    else:
185        return 1.0
186               
187def Ruland(RulCoff,wave,Q,Compton):
188    C = 2.9978e8
189    D = 1.5e-3
190    hmc = 0.024262734687
191    sinth2 = (Q*wave/(4.0*np.pi))**2
192    dlam = (wave**2)*Compton*Q/C
193    dlam_c = 2.0*hmc*sinth2-D*wave**2
194    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
195   
196def LorchWeight(Q):
197    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
198           
199def GetAsfMean(ElList,Sthl2):
200#   Calculate various scattering factor terms for PDF calcs
201#   ElList: element dictionary contains scattering factor coefficients, etc.
202#   Sthl2: numpy array of sin theta/lambda squared values
203#   returns: mean(f^2), mean(f)^2, mean(compton)
204    sumNoAtoms = 0.0
205    FF = np.zeros_like(Sthl2)
206    FF2 = np.zeros_like(Sthl2)
207    CF = np.zeros_like(Sthl2)
208    for El in ElList:
209        sumNoAtoms += ElList[El]['FormulaNo']
210    for El in ElList:
211        el = ElList[El]
212        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
213        cf = G2elem.ComptonFac(el,Sthl2)
214        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
215        FF2 += ff2*el['FormulaNo']/sumNoAtoms
216        CF += cf*el['FormulaNo']/sumNoAtoms
217    return FF2,FF**2,CF
218   
219def GetNumDensity(ElList,Vol):
220    sumNoAtoms = 0.0
221    for El in ElList:
222        sumNoAtoms += ElList[El]['FormulaNo']
223    return sumNoAtoms/Vol
224           
225def MultGetQ(Tth,MuT,Geometry,b=88.0,a=0.01):
226    NS = 500
227    Gama = np.linspace(0.,np.pi/2.,NS,False)[1:]
228    Cgama = np.cos(Gama)[:,np.newaxis]
229    Sgama = np.sin(Gama)[:,np.newaxis]
230    CSgama = 1.0/Sgama
231    Delt = Gama[1]-Gama[0]
232    emc = 7.94e-26
233    Navo = 6.023e23
234    Cth = npcosd(Tth/2.0)
235    CTth = npcosd(Tth)
236    Sth = npcosd(Tth/2.0)
237    STth = npsind(Tth)
238    CSth = 1.0/Sth
239    CSTth = 1.0/STth
240    SCth = 1.0/Cth
241    SCTth = 1.0/CTth
242    if 'Bragg' in Geometry:
243        G = 8.0*Delt*Navo*emc*Sth/((1.0-CTth**2)*(1.0-np.exp(-2.0*MuT*CSth)))
244        Y1 = np.pi
245        Y2 = np.pi/2.0
246        Y3 = 3.*np.pi/8. #3pi/4?
247        W = 1.0/(Sth+np.fabs(Sgama))
248        W += np.exp(-MuT*CSth)*(2.0*np.fabs(Sgama)*np.exp(-MuT*np.fabs(CSgama))-
249            (Sth+np.fabs(Sgama))*np.exp(-MuT*CSth))/(Sth**2-Sgama**2)
250        Fac0 = Sth**2*Sgama**2
251        X = Fac0+(Fac0+CTth)**2/2
252        Y = Cgama**2*Cth**2*(1.0-Fac0-CTth)
253        Z = Cgama**4*Cth**4/2.0
254        E = 2.0*(1.0-a)/(b*Cgama/Cth)
255        F1 = (2.0+b*(1.0+Sth*Sgama))/(b*Cth*Cgama) #trouble if < 1
256        F2 = (2.0+b*(1.0-Sth*Sgama))/(b*Cth*Cgama)
257        T1 = np.pi/np.sqrt(F1**2-1.0)
258        T2 = np.pi/np.sqrt(F2**2-1.0)
259        Y4 = T1+T2
260        Y5 = F1**2*T1+F2**2*T2-np.pi*(F1+F2)
261        Y6 = F1**4*T1+F2**4*T2-np.pi*(F1+F2)/2.0-np.pi*(F1**3+F2**3)
262        Y7 = (T2-T1)/(F1-F2)
263        YT = F2**2*T2-F1**2*T1
264        Y8 = Y1+YT/(F1-F2)
265        Y9 = Y2+(F2**4*T2-F1**4*T1)/(F1-F2)+Y1*((F1+F2)**2-F1*F2)
266        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
267       
268        Q = np.sum(W*M,axis=0)
269        return Q*G*NS/(NS-1.)
270#
271#      cos2th=2.0d*costh^2 - 1.0d
272#      G= delta * 8.0d * Na * emc * sinth/(1.0d + cos2th^2)/(1.0d - exp(-2.0d*mut*cscth))
273#      Y1=3.1415926d
274#      Y2=Y1*0.5d
275#      Y3=Y2*0.75d
276#      for i=1,num_steps-1 do begin
277#         cosgama=double(cos(gama[i]))
278#         singama=double(sin(gama[i]))
279#         cscgama=1.0d / singama
280#
281#         W=1.0d/(sinth+abs(singama))
282#         W=W+exp(-1.0*mut*cscth)*(2.0d*abs(singama)*exp(-1.0d*mut*abs(cscgama))- $
283#                 (sinth+abs(singama))*exp(-1.0d*mut*cscth))/(sinth^2-singama^2)
284#
285#         factor0=sinth^2*singama^2
286#         X=factor0+(factor0+cos2th)^2/2.0d
287#         Y=cosgama^2*(1.0d - factor0-cos2th)*costh^2
288#         Z=cosgama^4/2.0d*costh^4
289#         E=2.0d*(1.0-a)/b/cosgama/costh
290#
291#         F1=1.0d/b/cosgama*(2.0d + b*(1.0+sinth*singama))/costh
292#         F2=1.0d/b/cosgama*(2.0d + b*(1.0-sinth*singama))/costh
293#
294#         T1=3.14159/sqrt(F1^2-1.0d)
295#         T2=3.14159/sqrt(F2^2-1.0d)
296#         Y4=T1+T2
297#         Y5=F1^2*T1+F2^2*T2-3.14159*(F1+F2)
298#         Y6=F1^4*T1+F2^4*T2-3.14159*(F1+F2)/2.0-3.14159*(F1^3+F2^3)
299#         Y7=(T2-T1)/(F1-F2)
300#         Y8=Y1+(F2^2*T2-F1^2*T1)/(F1-F2)
301#         Y9=Y2+(F2^4*T2-F1^4*T1)/(F1-F2)+Y1*((F1+F2)^2-F1*F2)
302#         M=(a^2*(X*Y1+Y*Y2+Z*Y3)+a*E*(X*Y4+Y*Y5+Z*Y6)+E^2* $
303#                      (X*Y7+Y*Y8+Z*Y9))*cosgama
304#
305#         Q=Q+W*M
306#
307#      endfor
308#      Q=double(num_steps)/Double(num_steps-1)*Q*G
309#      end
310    elif 'Cylinder' in Geometry:
311        Q = 0.
312    elif 'Fixed' in Geometry:   #Dwiggens & Park, Acta Cryst. A27, 264 (1971) with corrections
313        EMA = np.exp(-MuT*(1.0-SCTth))
314        Fac1 = (1.-EMA)/(1.0-SCTth)
315        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*Fac1)
316        Fac0 = Cgama/(1-Sgama*SCTth)
317        Wp = Fac0*(Fac1-(EMA-np.exp(-MuT*(CSgama-SCTth)))/(CSgama-1.0))
318        Fac0 = Cgama/(1.0+Sgama*SCTth)
319        Wm = Fac0*(Fac1+(np.exp(-MuT*(1.0+CSgama))-1.0)/(CSgama+1.0))
320        X = (Sgama**2+CTth**2*(1.0-Sgama**2+Sgama**4))/2.0
321        Y = Sgama**3*Cgama*STth*CTth
322        Z = Cgama**2*(1.0+Sgama**2)*STth**2/2.0
323        Fac2 = 1.0/(b*Cgama*STth)
324        U = 2.0*(1.0-a)*Fac2
325        V = (2.0+b*(1.0-CTth*Sgama))*Fac2
326        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)
327        V = (2.0+b*(1.0+CTth*Sgama))*Fac2
328        Y = -Y
329        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)
330        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
331        return Q*G*NS/(NS-1.)
332    elif 'Tilting' in Geometry:
333        EMA = np.exp(-MuT*SCth)
334        G = 2.0*Delt*Navo*emc/((1.0+CTth**2)*EMA)
335#          Wplus = expmutsecth/(1.0d - singama*secth) + singama/mut/(1.0 -singama*secth)/(1.0-singama*secth)* $
336#                                                       (Exp(-1.0*mut*cscgama) - expmutsecth)
337#          Wminus = expmutsecth/(1.0d + singama*secth) - singama/mut/(1.0 +singama*secth)/(1.0+singama*secth)* $
338#                                                        expmutsecth*(1.0d - expmutsecth*Exp(-1.0*mut*cscgama))
339        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)/(1.0-Sgama*SCth)*(np.exp(-MuT*CSgama)-EMA)
340#        Wp = EMA/(1.0-Sgama*SCth)+Sgama/MuT/(1.0-Sgama*SCth)**2*(np.exp(-MuT*CSgama)-EMA)
341        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)/(1.0+Sgama*SCth)*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
342#        Wm = EMA/(1.0+Sgama*SCth)-Sgama/MuT/(1.0+Sgama*SCth)**2*EMA*(1.0-EMA*np.exp(-MuT*CSgama))
343        X = 0.5*(Cth**2*(Cth**2*Sgama**4-4.0*Sth**2*Cgama**2)+1.0)
344        Y = Cgama**2*(1.0+Cgama**2)*Cth**2*Sth**2
345        Z = 0.5*Cgama**4*Sth**4
346#          X = 0.5*(costh*costh*(costh*costh*singama*singama*singama*singama - $
347#                           4.0*sinth*sinth*cosgama*cosgama) +1.0d)
348#
349#          Y = cosgama*cosgama*(1.0 + cosgama*cosgama)*costh*costh*sinth*sinth
350#
351#          Z= 0.5 *cosgama*cosgama*cosgama*cosgama* (sinth^4)
352#
353        AlP = 2.0+b*(1.0-Cth*Sgama)
354        AlM = 2.0+b*(1.0+Cth*Sgama)
355#          alphaplus = 2.0 + b*(1.0 - costh*singama)
356#          alphaminus = 2.0 + b*(1.0 + costh*singama)
357        BeP = np.sqrt(np.fabs(AlP**2-(b*Cgama*Sth)**2))
358        BeM = np.sqrt(np.fabs(AlM**2-(b*Cgama*Sth)**2))
359#          betaplus = Sqrt(Abs(alphaplus*alphaplus - b*b*cosgama*cosgama*sinth*sinth))
360#          betaminus = Sqrt(Abs(alphaminus*alphaminus - b*b*cosgama*cosgama*sinth*sinth))
361        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)* \
362            (4.0*X/AlP/BeP+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlP/BeP-1.0)+
363            2.0/b**4*AlP/BeP*AlP**2-2.0/b**4*AlP**2-Cgama**2/b/b*Sth*2))
364#          Mplus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
365#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaplus)* $
366#                   (4.0*x/alphaplus/betaplus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaplus/betaplus -1.0) + $
367#                   2.0/(b^4)*alphaplus/betaplus*alphaplus*alphaplus - 2.0/(b^4)*alphaplus*alphaplus - $
368#                   cosgama*cosgama/b/b*sinth*sinth))
369        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)* \
370            (4.0*X/AlM/BeM+(4.0*(1.0+Cgama**2)/b/b*Cth**2)*(AlM/BeM-1.0)+
371            2.0/b**4*AlM/BeM*AlM**2-2.0/b**4*AlM**2-Cgama**2/b/b*Sth*2))
372#          Mminus = cosgama*(!DPI * a * a * (2.0*x + y + 0.75*z) + $
373#                   (2.0*!DPI*(1.0 - a)) *(1.0 - a + a*alphaminus)* $
374#                   (4.0*x/alphaminus/betaminus + (4.0*(1.0+cosgama*cosgama)/b/b*costh*costh)*(alphaminus/betaminus -1.0) + $
375#                   2.0/(b^4)*alphaminus/betaminus*alphaminus*alphaminus - 2.0/(b^4)*alphaminus*alphaminus - $
376#                   cosgama*cosgama/b/b*sinth*sinth))
377        Q = np.sum(Wp*Mp+Wm*Mm,axis=0)
378        return Q*G*NS/(NS-1.)
379#       expmutsecth = Exp(-1.0*mut*secth)
380#       G= delta * 2.0 * Na * emc /(1.0+costth^2)/expmutsecth
381#       for i=1, num_steps-1 do begin
382#          cosgama=double(cos(gama[i]))
383#          singama=double(sin(gama[i]))
384#          cscgama=1.0d/singama
385#
386#
387#     ; print, "W", min(wplus), max(wplus), min(wminus), max(wminus)
388#
389#
390#
391#
392#    ;               print, a,b
393#  ; print, "M", min(mplus), max(mplus), min(mminus), max(mminus)
394#          Q=Q+ Wplus*Mplus + Wminus*Mminus
395#      endfor
396#      Q=double(num_steps)/double(num_steps-1)*Q*G
397#   ;   print, min(q), max(q)
398#      end
399
400def MultiScattering(Geometry,ElList,Tth):
401    BN = BD = 0.0
402    Amu = 0.0
403    for El in ElList:
404        el = ElList[El]
405        BN += el['Z']*el['FormulaNo']
406        BD += el['FormulaNo']
407        Amu += el['FormulaNo']*el['mu']
408       
409#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
410
411np.seterr(divide='ignore')
412
413# Normal distribution
414
415# loc = mu, scale = std
416_norm_pdf_C = 1./math.sqrt(2*math.pi)
417class norm_gen(st.rv_continuous):
418       
419    def pdf(self,x,*args,**kwds):
420        loc,scale=kwds['loc'],kwds['scale']
421        x = (x-loc)/scale
422        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
423       
424norm = norm_gen(name='norm',longname='A normal',extradoc="""
425
426Normal distribution
427
428The location (loc) keyword specifies the mean.
429The scale (scale) keyword specifies the standard deviation.
430
431normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
432""")
433
434## Cauchy
435
436# median = loc
437
438class cauchy_gen(st.rv_continuous):
439
440    def pdf(self,x,*args,**kwds):
441        loc,scale=kwds['loc'],kwds['scale']
442        x = (x-loc)/scale
443        return 1.0/np.pi/(1.0+x*x) / scale
444       
445cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
446
447Cauchy distribution
448
449cauchy.pdf(x) = 1/(pi*(1+x**2))
450
451This is the t distribution with one degree of freedom.
452""")
453   
454   
455#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
456
457
458class fcjde_gen(st.rv_continuous):
459    """
460    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
461    Ref: J. Appl. Cryst. (1994) 27, 892-900.
462    Parameters
463    -----------------------------------------
464    x: array -1 to 1
465    t: 2-theta position of peak
466    s: sum(S/L,H/L); S: sample height, H: detector opening,
467        L: sample to detector opening distance
468    dx: 2-theta step size in deg
469    Result for fcj.pdf
470    -----------------------------------------
471    T = x*dx+t
472    s = S/L+H/L
473    if x < 0:
474        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
475    if x >= 0:
476        fcj.pdf = 0   
477    """
478    def _pdf(self,x,t,s,dx):
479        T = dx*x+t
480        ax2 = abs(npcosd(T))
481        ax = ax2**2
482        bx = npcosd(t)**2
483        bx = np.where(ax>bx,bx,ax)
484        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
485        fx = np.where(fx > 0.,fx,0.0)
486        return fx
487             
488    def pdf(self,x,*args,**kwds):
489        loc=kwds['loc']
490        return self._pdf(x-loc,*args)
491       
492fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
493               
494def getBackground(pfx,parmDict,bakType,xdata):
495    yb = np.zeros_like(xdata)
496    if bakType == 'chebyschev':
497        iBak = 0
498        while True:
499            key = pfx+'Back:'+str(iBak)
500            try:
501                yb += parmDict[key]*(xdata-xdata[0])**iBak
502                iBak += 1
503            except KeyError:
504                break
505    return yb
506
507def getWidths(pos,sig,gam,shl):
508    widths = [np.sqrt(sig)/100.,gam/200.]
509    fwhm = 2.355*widths[0]+2.*widths[1]
510    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
511    fmax = 15.0*fwhm
512    if pos > 90:
513        fmin,fmax = [fmax,fmin]         
514    return widths,fmin,fmax
515               
516def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
517    DX = xdata[1]-xdata[0]
518    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
519    x = np.linspace(pos-fmin,pos+fmin,256)
520    dx = x[1]-x[0]
521    Norm = norm.pdf(x,loc=pos,scale=widths[0])
522    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
523    arg = [pos,shl/57.2958,dx,]
524    FCJ = fcjde.pdf(x,*arg,loc=pos)
525    if len(np.nonzero(FCJ)[0])>5:
526        z = np.column_stack([Norm,Cauchy,FCJ]).T
527        Z = fft(z)
528        Df = ifft(Z.prod(axis=0)).real
529    else:
530        z = np.column_stack([Norm,Cauchy]).T
531        Z = fft(z)
532        Df = fftshift(ifft(Z.prod(axis=0))).real
533    Df /= np.sum(Df)
534    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
535    return intens*Df(xdata)*DX/dx
536
537#use old fortran routine
538def getFCJVoigt3(pos,intens,sig,gam,shl,xdata):
539   
540    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
541    Df /= np.sum(Df)
542    return intens*Df
543
544def getPeakProfile(parmDict,xdata,varyList,bakType):
545   
546    yb = getBackground('',parmDict,bakType,xdata)
547    yc = np.zeros_like(yb)
548    dx = xdata[1]-xdata[0]
549    U = parmDict['U']
550    V = parmDict['V']
551    W = parmDict['W']
552    X = parmDict['X']
553    Y = parmDict['Y']
554    shl = max(parmDict['SH/L'],0.0005)
555    Ka2 = False
556    if 'Lam1' in parmDict.keys():
557        Ka2 = True
558        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
559        kRatio = parmDict['I(L2)/I(L1)']
560    iPeak = 0
561    while True:
562        try:
563            pos = parmDict['pos'+str(iPeak)]
564            intens = parmDict['int'+str(iPeak)]
565            sigName = 'sig'+str(iPeak)
566            if sigName in varyList:
567                sig = parmDict[sigName]
568            else:
569                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
570            sig = max(sig,0.001)          #avoid neg sigma
571            gamName = 'gam'+str(iPeak)
572            if gamName in varyList:
573                gam = parmDict[gamName]
574            else:
575                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
576            gam = max(gam,0.001)             #avoid neg gamma
577            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
578            iBeg = np.searchsorted(xdata,pos-fmin)
579            lenX = len(xdata)
580            if not iBeg:
581                iFin = np.searchsorted(xdata,pos+fmin)
582            elif iBeg == lenX:
583                iFin = iBeg
584            else:
585                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
586            if not iBeg+iFin:       #peak below low limit
587                iPeak += 1
588                continue
589            elif not iBeg-iFin:     #peak above high limit
590                return yb+yc
591            yc[iBeg:iFin] += getFCJVoigt3(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
592            if Ka2:
593                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
594                kdelt = int((pos2-pos)/dx)               
595                iBeg = min(lenX,iBeg+kdelt)
596                iFin = min(lenX,iFin+kdelt)
597                if iBeg-iFin:
598                    yc[iBeg:iFin] += getFCJVoigt3(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
599            iPeak += 1
600        except KeyError:        #no more peaks to process
601            return yb+yc
602
603def Dict2Values(parmdict, varylist):
604    '''Use before call to leastsq to setup list of values for the parameters
605    in parmdict, as selected by key in varylist'''
606    return [parmdict[key] for key in varylist] 
607   
608def Values2Dict(parmdict, varylist, values):
609    ''' Use after call to leastsq to update the parameter dictionary with
610    values corresponding to keys in varylist'''
611    parmdict.update(zip(varylist,values))
612   
613def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
614   
615    def SetBackgroundParms(Background):
616        bakType,bakFlag = Background[:2]
617        backVals = Background[3:]
618        backNames = ['Back:'+str(i) for i in range(len(backVals))]
619        if bakFlag: #returns backNames as varyList = backNames
620            return bakType,dict(zip(backNames,backVals)),backNames
621        else:       #no background varied; varyList = []
622            return bakType,dict(zip(backNames,backVals)),[]
623       
624    def GetBackgroundParms(parmList,Background):
625        iBak = 0
626        while True:
627            try:
628                bakName = 'Back:'+str(iBak)
629                Background[iBak+3] = parmList[bakName]
630                iBak += 1
631            except KeyError:
632                break
633               
634    def BackgroundPrint(Background,sigDict):
635        if Background[1]:
636            print 'Background coefficients for',Background[0],'function'
637            ptfmt = "%12.5f"
638            ptstr =  'values:'
639            sigstr = 'esds  :'
640            for i,back in enumerate(Background[3:]):
641                ptstr += ptfmt % (back)
642                sigstr += ptfmt % (sigDict['Back:'+str(i)])
643            print ptstr
644            print sigstr
645        else:
646            print 'Background not refined'
647           
648    def SetInstParms(Inst):
649        insVals,insFlags,insNames = Inst[1:4]
650        dataType = insVals[0]
651        insVary = []
652        for i,flag in enumerate(insFlags):
653            if flag:
654                insVary.append(insNames[i])
655        instDict = dict(zip(insNames,insVals))
656        instDict['X'] = max(instDict['X'],0.01)
657        instDict['Y'] = max(instDict['Y'],0.01)
658        instDict['SH/L'] = max(instDict['SH/L'],0.0001)
659        return dataType,instDict,insVary
660       
661    def GetInstParms(parmDict,Inst,varyList,Peaks):
662        instNames = Inst[3]
663        for i,name in enumerate(instNames):
664            Inst[1][i] = parmDict[name]
665        iPeak = 0
666        while True:
667            try:
668                sigName = 'sig'+str(iPeak)
669                pos = parmDict['pos'+str(iPeak)]
670                if sigName not in varyList:
671                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
672                gamName = 'gam'+str(iPeak)
673                if gamName not in varyList:
674                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
675                iPeak += 1
676            except KeyError:
677                break
678       
679    def InstPrint(Inst,sigDict):
680        print 'Instrument Parameters:'
681        ptfmt = "%12.6f"
682        ptlbls = 'names :'
683        ptstr =  'values:'
684        sigstr = 'esds  :'
685        instNames = Inst[3][1:]
686        for i,name in enumerate(instNames):
687            ptlbls += "%s" % (name.center(12))
688            ptstr += ptfmt % (Inst[1][i+1])
689            if name in sigDict:
690                sigstr += ptfmt % (sigDict[name])
691            else:
692                sigstr += 12*' '
693        print ptlbls
694        print ptstr
695        print sigstr
696
697    def SetPeaksParms(Peaks):
698        peakNames = []
699        peakVary = []
700        peakVals = []
701        names = ['pos','int','sig','gam']
702        for i,peak in enumerate(Peaks):
703            for j in range(4):
704                peakVals.append(peak[2*j])
705                parName = names[j]+str(i)
706                peakNames.append(parName)
707                if peak[2*j+1]:
708                    peakVary.append(parName)
709        return dict(zip(peakNames,peakVals)),peakVary
710               
711    def GetPeaksParms(parmDict,Peaks,varyList):
712        names = ['pos','int','sig','gam']
713        for i,peak in enumerate(Peaks):
714            for j in range(4):
715                pos = parmDict['pos'+str(i)]
716                parName = names[j]+str(i)
717                if parName in varyList:
718                    peak[2*j] = parmDict[parName]
719                elif 'sig' in parName:
720                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
721                elif 'gam' in parName:
722                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
723                       
724    def PeaksPrint(parmDict,sigDict,varyList):
725        print 'Peak coefficients:'
726        names = ['pos','int','sig','gam']
727        head = 15*' '
728        for name in names:
729            head += name.center(12)+'esd'.center(12)
730        print head
731        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
732        for i,peak in enumerate(Peaks):
733            ptstr =  ':'
734            for j in range(4):
735                name = names[j]
736                parName = name+str(i)
737                ptstr += ptfmt[name] % (parmDict[parName])
738                if parName in varyList:
739#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
740                    ptstr += ptfmt[name] % (sigDict[parName])
741                else:
742#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
743                    ptstr += 12*' '
744            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
745           
746    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
747        parmdict.update(zip(varylist,values))
748        M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
749        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
750        if dlg:
751            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
752            if not GoOn:
753                return -M           #abort!!
754        return M
755       
756    Ftol = 0.0001
757    if oneCycle:
758        Ftol = 1.0
759    x,y,w,yc,yb,yd = data               #these are numpy arrays!
760    xBeg = np.searchsorted(x,Limits[0])
761    xFin = np.searchsorted(x,Limits[1])
762    bakType,bakDict,bakVary = SetBackgroundParms(Background)
763    dataType,insDict,insVary = SetInstParms(Inst)
764    peakDict,peakVary = SetPeaksParms(Peaks)
765    parmDict = {}
766    parmDict.update(bakDict)
767    parmDict.update(insDict)
768    parmDict.update(peakDict)
769    varyList = bakVary+insVary+peakVary
770    while True:
771        begin = time.time()
772        values =  np.array(Dict2Values(parmDict, varyList))
773        if FitPgm == 'LSQ':
774            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
775            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
776            screenSize = wx.ClientDisplayRect()
777            Size = dlg.GetSize()
778            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
779            try:
780                result = so.leastsq(errPeakProfile,values,full_output=True,epsfcn=1.e-8,ftol=Ftol,
781                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
782            finally:
783                dlg.Destroy()
784            runtime = time.time()-begin   
785            chisq = np.sum(result[2]['fvec']**2)
786            ncyc = int(result[2]['nfev']/len(varyList))
787            Values2Dict(parmDict, varyList, result[0])
788            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
789            GOF = chisq/(xFin-xBeg-len(varyList))
790            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
791            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
792            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
793            try:
794                sig = np.sqrt(np.diag(result[1])*GOF)
795                if np.any(np.isnan(sig)):
796                    print '*** Least squares aborted - some invalid esds possible ***'
797                break                   #refinement succeeded - finish up!
798            except ValueError:          #result[1] is None on singular matrix
799                print '**** Refinement failed - singular matrix ****'
800                Ipvt = result[2]['ipvt']
801                for i,ipvt in enumerate(Ipvt):
802                    if not np.sum(result[2]['fjac'],axis=1)[i]:
803                        print 'Removing parameter: ',varyList[ipvt-1]
804                        del(varyList[ipvt-1])
805                        break
806        elif FitPgm == 'BFGS':
807            print 'Other program here'
808            return
809       
810    sigDict = dict(zip(varyList,sig))
811    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
812    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
813    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
814    GetBackgroundParms(parmDict,Background)
815    BackgroundPrint(Background,sigDict)
816    GetInstParms(parmDict,Inst,varyList,Peaks)
817    InstPrint(Inst,sigDict)
818    GetPeaksParms(parmDict,Peaks,varyList)   
819    PeaksPrint(parmDict,sigDict,varyList)
820   
821def CalcPDF(data,inst,xydata):
822    auxPlot = []
823    import copy
824    import scipy.fftpack as ft
825    #subtract backgrounds - if any
826    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
827    if data['Sample Bkg.']['Name']:
828        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
829            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
830    if data['Container']['Name']:
831        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
832        if data['Container Bkg.']['Name']:
833            xycontainer += (xydata['Container Bkg.'][1][1]+
834                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
835        xydata['IofQ'][1][1] += xycontainer
836    #get element data & absorption coeff.
837    ElList = data['ElList']
838    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
839    #Apply angle dependent corrections
840    Tth = xydata['Sample'][1][0]
841    dt = (Tth[1]-Tth[0])
842    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
843    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
844    if data['DetType'] == 'Image plate':
845        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
846    XY = xydata['IofQ'][1]   
847    #convert to Q
848    hc = 12.397639
849    if 'Lam' in inst:
850        wave = inst['Lam']
851    else:
852        wave = inst['Lam1']
853    keV = hc/wave
854    minQ = npT2q(Tth[0],wave)
855    maxQ = npT2q(Tth[-1],wave)   
856    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
857    dq = Qpoints[1]-Qpoints[0]
858    XY[0] = npT2q(XY[0],wave)   
859#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
860    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
861    Qdata = T(Qpoints)
862   
863    qLimits = data['QScaleLim']
864    minQ = np.searchsorted(Qpoints,qLimits[0])
865    maxQ = np.searchsorted(Qpoints,qLimits[1])
866    newdata = []
867    xydata['IofQ'][1][0] = Qpoints
868    xydata['IofQ'][1][1] = Qdata
869    for item in xydata['IofQ'][1]:
870        newdata.append(item[:maxQ])
871    xydata['IofQ'][1] = newdata
872   
873
874    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
875    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
876    Q = xydata['SofQ'][1][0]
877    ruland = Ruland(data['Ruland'],wave,Q,CF)
878#    auxPlot.append([Q,ruland,'Ruland'])     
879    CF *= ruland
880#    auxPlot.append([Q,CF,'CF-Corr'])
881    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
882    xydata['SofQ'][1][1] *= scale
883    xydata['SofQ'][1][1] -= CF
884    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
885    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
886    xydata['SofQ'][1][1] *= scale
887   
888    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
889    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
890    if data['Lorch']:
891        xydata['FofQ'][1][1] *= LorchWeight(Q)
892   
893    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
894    nR = len(xydata['GofR'][1][1])
895    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
896    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
897   
898       
899    return auxPlot
900       
901#testing data
902NeedTestData = True
903def TestData():
904#    global NeedTestData
905    NeedTestData = False
906    global bakType
907    bakType = 'chebyschev'
908    global xdata
909    xdata = np.linspace(4.0,40.0,36000)
910    global parmDict0
911    parmDict0 = {
912        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
913        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
914        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
915        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
916        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
917        'Back0':5.384,'Back1':-0.015,'Back2':.004,
918        }
919    global parmDict1
920    parmDict1 = {
921        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
922        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
923        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
924        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
925        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
926        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
927        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
928        'Back0':36.897,'Back1':-0.508,'Back2':.006,
929        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
930        }
931    global varyList
932    varyList = []
933
934def test0():
935    if NeedTestData: TestData()
936    msg = 'test '
937    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
938    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
939    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
940    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
941    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
942    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
943   
944def test1():
945    if NeedTestData: TestData()
946    time0 = time.time()
947    for i in range(100):
948        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
949    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
950   
951   
952
953if __name__ == '__main__':
954    import GSASIItestplot as plot
955    global plotter
956    plotter = plot.PlotNotebook()
957    test0()
958    print "OK"
959    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.