source: trunk/GSASIIpwd.py @ 342

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

major modifications to get 1st "working" version of Refine
Add GSASIImapvars.py
fix cell refinement in indexing window
single cycle for peak refinement

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