source: trunk/GSASIIpwd.py @ 345

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

put new banner in GSASII.BAT
fix for vertical tilt & offset images
speedup of peak calcs
allow change from lam to lam1 & lam2
Pawley refinement fixes

File size: 36.7 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# Normal distribution
413
414# loc = mu, scale = std
415_norm_pdf_C = 1./math.sqrt(2*math.pi)
416class norm_gen(st.rv_continuous):
417       
418    def pdf(self,x,*args,**kwds):
419        loc,scale=kwds['loc'],kwds['scale']
420        x = (x-loc)/scale
421        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
422       
423norm = norm_gen(name='norm',longname='A normal',extradoc="""
424
425Normal distribution
426
427The location (loc) keyword specifies the mean.
428The scale (scale) keyword specifies the standard deviation.
429
430normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
431""")
432
433## Cauchy
434
435# median = loc
436
437class cauchy_gen(st.rv_continuous):
438
439    def pdf(self,x,*args,**kwds):
440        loc,scale=kwds['loc'],kwds['scale']
441        x = (x-loc)/scale
442        return 1.0/np.pi/(1.0+x*x) / scale
443       
444cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
445
446Cauchy distribution
447
448cauchy.pdf(x) = 1/(pi*(1+x**2))
449
450This is the t distribution with one degree of freedom.
451""")
452   
453   
454#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
455
456
457class fcjde_gen(st.rv_continuous):
458    """
459    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
460    Ref: J. Appl. Cryst. (1994) 27, 892-900.
461    Parameters
462    -----------------------------------------
463    x: array -1 to 1
464    t: 2-theta position of peak
465    s: sum(S/L,H/L); S: sample height, H: detector opening,
466        L: sample to detector opening distance
467    dx: 2-theta step size in deg
468    Result for fcj.pdf
469    -----------------------------------------
470    T = x*dx+t
471    s = S/L+H/L
472    if x < 0:
473        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
474    if x >= 0:
475        fcj.pdf = 0   
476    """
477    def _pdf(self,x,t,s,dx):
478        T = dx*x+t
479        ax2 = abs(npcosd(T))
480        ax = ax2**2
481        bx = npcosd(t)**2
482        bx = np.where(ax>bx,bx,ax)
483        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
484        fx = np.where(fx > 0.,fx,0.0)
485        return fx
486             
487    def pdf(self,x,*args,**kwds):
488        loc=kwds['loc']
489        return self._pdf(x-loc,*args)
490       
491fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
492               
493def getBackground(pfx,parmDict,bakType,xdata):
494    yb = np.zeros_like(xdata)
495    if bakType == 'chebyschev':
496        iBak = 0
497        while True:
498            key = pfx+'Back:'+str(iBak)
499            try:
500                yb += parmDict[key]*(xdata-xdata[0])**iBak
501                iBak += 1
502            except KeyError:
503                break
504    return yb
505
506def getWidths(pos,sig,gam,shl):
507    widths = [np.sqrt(sig)/100.,gam/200.]
508    fwhm = 2.355*widths[0]+2.*widths[1]
509    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
510    fmax = 15.0*fwhm
511    if pos > 90:
512        fmin,fmax = [fmax,fmin]         
513    return widths,fmin,fmax
514               
515def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
516    DX = xdata[1]-xdata[0]
517    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
518    x = np.linspace(pos-fmin,pos+fmin,256)
519    dx = x[1]-x[0]
520    Norm = norm.pdf(x,loc=pos,scale=widths[0])
521    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
522    arg = [pos,shl/57.2958,dx,]
523    FCJ = fcjde.pdf(x,*arg,loc=pos)
524    if len(np.nonzero(FCJ)[0])>5:
525        z = np.column_stack([Norm,Cauchy,FCJ]).T
526        Z = fft(z)
527        Df = ifft(Z.prod(axis=0)).real
528    else:
529        z = np.column_stack([Norm,Cauchy]).T
530        Z = fft(z)
531        Df = fftshift(ifft(Z.prod(axis=0))).real
532    Df /= np.sum(Df)
533    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
534    return intens*Df(xdata)*DX/dx
535
536def getPeakProfile(parmDict,xdata,varyList,bakType):
537   
538    yb = getBackground('',parmDict,bakType,xdata)
539    yc = np.zeros_like(yb)
540    dx = xdata[1]-xdata[0]
541    U = parmDict['U']
542    V = parmDict['V']
543    W = parmDict['W']
544    X = parmDict['X']
545    Y = parmDict['Y']
546    shl = max(parmDict['SH/L'],0.0005)
547    Ka2 = False
548    if 'Lam1' in parmDict.keys():
549        Ka2 = True
550        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
551        kRatio = parmDict['I(L2)/I(L1)']
552    iPeak = 0
553    while True:
554        try:
555            pos = parmDict['pos'+str(iPeak)]
556            intens = parmDict['int'+str(iPeak)]
557            sigName = 'sig'+str(iPeak)
558            if sigName in varyList:
559                sig = parmDict[sigName]
560            else:
561                sig = U*tand(pos/2.0)**2+V*tand(pos/2.0)+W
562            sig = max(sig,0.001)          #avoid neg sigma
563            gamName = 'gam'+str(iPeak)
564            if gamName in varyList:
565                gam = parmDict[gamName]
566            else:
567                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
568            gam = max(gam,0.001)             #avoid neg gamma
569            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
570            iBeg = np.searchsorted(xdata,pos-fmin)
571            lenX = len(xdata)
572            if not iBeg:
573                iFin = np.searchsorted(xdata,pos+fmin)
574            elif iBeg == lenX:
575                iFin = iBeg
576            else:
577                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
578            if not iBeg+iFin:       #peak below low limit
579                iPeak += 1
580                continue
581            elif not iBeg-iFin:     #peak above high limit
582                return yb+yc
583            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
584            if Ka2:
585                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
586                kdelt = int((pos2-pos)/dx)               
587                iBeg = min(lenX,iBeg+kdelt)
588                iFin = min(lenX,iFin+kdelt)
589                if iBeg-iFin:
590                    yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
591            iPeak += 1
592        except KeyError:        #no more peaks to process
593            return yb+yc
594
595def Dict2Values(parmdict, varylist):
596    '''Use before call to leastsq to setup list of values for the parameters
597    in parmdict, as selected by key in varylist'''
598    return [parmdict[key] for key in varylist] 
599   
600def Values2Dict(parmdict, varylist, values):
601    ''' Use after call to leastsq to update the parameter dictionary with
602    values corresponding to keys in varylist'''
603    parmdict.update(zip(varylist,values))
604   
605def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,data,oneCycle=False):
606   
607    def SetBackgroundParms(Background):
608        bakType,bakFlag = Background[:2]
609        backVals = Background[3:]
610        backNames = ['Back:'+str(i) for i in range(len(backVals))]
611        if bakFlag: #returns backNames as varyList = backNames
612            return bakType,dict(zip(backNames,backVals)),backNames
613        else:       #no background varied; varyList = []
614            return bakType,dict(zip(backNames,backVals)),[]
615       
616    def GetBackgroundParms(parmList,Background):
617        iBak = 0
618        while True:
619            try:
620                bakName = 'Back:'+str(iBak)
621                Background[iBak+3] = parmList[bakName]
622                iBak += 1
623            except KeyError:
624                break
625               
626    def BackgroundPrint(Background,sigDict):
627        if Background[1]:
628            print 'Background coefficients for',Background[0],'function'
629            ptfmt = "%12.5f"
630            ptstr =  'values:'
631            sigstr = 'esds  :'
632            for i,back in enumerate(Background[3:]):
633                ptstr += ptfmt % (back)
634                sigstr += ptfmt % (sigDict['Back:'+str(i)])
635            print ptstr
636            print sigstr
637        else:
638            print 'Background not refined'
639           
640    def SetInstParms(Inst):
641        insVals,insFlags,insNames = Inst[1:4]
642        dataType = insVals[0]
643        insVary = []
644        for i,flag in enumerate(insFlags):
645            if flag:
646                insVary.append(insNames[i])
647        instDict = dict(zip(insNames,insVals))
648        instDict['X'] = max(instDict['X'],0.01)
649        instDict['Y'] = max(instDict['Y'],0.01)
650        instDict['SH/L'] = max(instDict['SH/L'],0.0001)
651        return dataType,instDict,insVary
652       
653    def GetInstParms(parmDict,Inst,varyList,Peaks):
654        instNames = Inst[3]
655        for i,name in enumerate(instNames):
656            Inst[1][i] = parmDict[name]
657        iPeak = 0
658        while True:
659            try:
660                sigName = 'sig'+str(iPeak)
661                pos = parmDict['pos'+str(iPeak)]
662                if sigName not in varyList:
663                    parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
664                gamName = 'gam'+str(iPeak)
665                if gamName not in varyList:
666                    parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
667                iPeak += 1
668            except KeyError:
669                break
670       
671    def InstPrint(Inst,sigDict):
672        print 'Instrument Parameters:'
673        ptfmt = "%12.6f"
674        ptlbls = 'names :'
675        ptstr =  'values:'
676        sigstr = 'esds  :'
677        instNames = Inst[3][1:]
678        for i,name in enumerate(instNames):
679            ptlbls += "%s" % (name.center(12))
680            ptstr += ptfmt % (Inst[1][i+1])
681            if name in sigDict:
682                sigstr += ptfmt % (sigDict[name])
683            else:
684                sigstr += 12*' '
685        print ptlbls
686        print ptstr
687        print sigstr
688
689    def SetPeaksParms(Peaks):
690        peakNames = []
691        peakVary = []
692        peakVals = []
693        names = ['pos','int','sig','gam']
694        for i,peak in enumerate(Peaks):
695            for j in range(4):
696                peakVals.append(peak[2*j])
697                parName = names[j]+str(i)
698                peakNames.append(parName)
699                if peak[2*j+1]:
700                    peakVary.append(parName)
701        return dict(zip(peakNames,peakVals)),peakVary
702               
703    def GetPeaksParms(parmDict,Peaks,varyList):
704        names = ['pos','int','sig','gam']
705        for i,peak in enumerate(Peaks):
706            for j in range(4):
707                pos = parmDict['pos'+str(i)]
708                parName = names[j]+str(i)
709                if parName in varyList:
710                    peak[2*j] = parmDict[parName]
711                elif 'sig' in parName:
712                    peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
713                elif 'gam' in parName:
714                    peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
715                       
716    def PeaksPrint(parmDict,sigDict,varyList):
717        print 'Peak coefficients:'
718        names = ['pos','int','sig','gam']
719        head = 15*' '
720        for name in names:
721            head += name.center(12)+'esd'.center(12)
722        print head
723        ptfmt = {'pos':"%12.5f",'int':"%12.1f",'sig':"%12.3f",'gam':"%12.3f"}
724        for i,peak in enumerate(Peaks):
725            ptstr =  ':'
726            for j in range(4):
727                name = names[j]
728                parName = name+str(i)
729                ptstr += ptfmt[name] % (parmDict[parName])
730                if parName in varyList:
731#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
732                    ptstr += ptfmt[name] % (sigDict[parName])
733                else:
734#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
735                    ptstr += 12*' '
736            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
737           
738    def errPeakProfile(values, xdata, ydata, weights, parmdict, varylist,bakType,dlg):       
739        parmdict.update(zip(varylist,values))
740        M = np.sqrt(weights)*(ydata-getPeakProfile(parmdict,xdata,varylist,bakType))
741        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
742        if dlg:
743            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
744            if not GoOn:
745                return -M           #abort!!
746        return M
747       
748    Ftol = 0.0001
749    if oneCycle:
750        Ftol = 1.0
751    x,y,w,yc,yb,yd = data               #these are numpy arrays!
752    xBeg = np.searchsorted(x,Limits[0])
753    xFin = np.searchsorted(x,Limits[1])
754    bakType,bakDict,bakVary = SetBackgroundParms(Background)
755    dataType,insDict,insVary = SetInstParms(Inst)
756    peakDict,peakVary = SetPeaksParms(Peaks)
757    parmDict = {}
758    parmDict.update(bakDict)
759    parmDict.update(insDict)
760    parmDict.update(peakDict)
761    varyList = bakVary+insVary+peakVary
762    while True:
763        begin = time.time()
764        values =  np.array(Dict2Values(parmDict, varyList))
765        if FitPgm == 'LSQ':
766            dlg = wx.ProgressDialog('Residual','Peak fit Rwp = ',101.0, 
767            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
768            screenSize = wx.ClientDisplayRect()
769            Size = dlg.GetSize()
770            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
771            try:
772                result = so.leastsq(errPeakProfile,values,full_output=True,             #ftol=Ftol,
773                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
774            finally:
775                dlg.Destroy()
776            runtime = time.time()-begin   
777            chisq = np.sum(result[2]['fvec']**2)
778            ncyc = int(result[2]['nfev']/len(varyList))
779            Values2Dict(parmDict, varyList, result[0])
780            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
781            GOF = chisq/(xFin-xBeg-len(varyList))
782            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
783            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
784            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
785            try:
786                sig = np.sqrt(np.diag(result[1])*GOF)
787                if np.any(np.isnan(sig)):
788                    print '*** Least squares aborted - some invalid esds possible ***'
789                break                   #refinement succeeded - finish up!
790            except ValueError:          #result[1] is None on singular matrix
791                print '**** Refinement failed - singular matrix ****'
792                Ipvt = result[2]['ipvt']
793                for i,ipvt in enumerate(Ipvt):
794                    if not np.sum(result[2]['fjac'],axis=1)[i]:
795                        print 'Removing parameter: ',varyList[ipvt-1]
796                        del(varyList[ipvt-1])
797                        break
798        elif FitPgm == 'BFGS':
799            print 'Other program here'
800            return
801       
802    sigDict = dict(zip(varyList,sig))
803    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
804    yc[xBeg:xFin] = getPeakProfile(parmDict,x[xBeg:xFin],varyList,bakType)
805    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
806    GetBackgroundParms(parmDict,Background)
807    BackgroundPrint(Background,sigDict)
808    GetInstParms(parmDict,Inst,varyList,Peaks)
809    InstPrint(Inst,sigDict)
810    GetPeaksParms(parmDict,Peaks,varyList)   
811    PeaksPrint(parmDict,sigDict,varyList)
812   
813def CalcPDF(data,inst,xydata):
814    auxPlot = []
815    import copy
816    import scipy.fftpack as ft
817    #subtract backgrounds - if any
818    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
819    if data['Sample Bkg.']['Name']:
820        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
821            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
822    if data['Container']['Name']:
823        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
824        if data['Container Bkg.']['Name']:
825            xycontainer += (xydata['Container Bkg.'][1][1]+
826                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
827        xydata['IofQ'][1][1] += xycontainer
828    #get element data & absorption coeff.
829    ElList = data['ElList']
830    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
831    #Apply angle dependent corrections
832    Tth = xydata['Sample'][1][0]
833    dt = (Tth[1]-Tth[0])
834    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],Abs,data['Diam'],Tth)
835    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'],Tth,Azm=inst['Azimuth'])
836    if data['DetType'] == 'Image plate':
837        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
838    XY = xydata['IofQ'][1]   
839    #convert to Q
840    hc = 12.397639
841    if 'Lam' in inst:
842        wave = inst['Lam']
843    else:
844        wave = inst['Lam1']
845    keV = hc/wave
846    minQ = npT2q(Tth[0],wave)
847    maxQ = npT2q(Tth[-1],wave)   
848    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
849    dq = Qpoints[1]-Qpoints[0]
850    XY[0] = npT2q(XY[0],wave)   
851#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
852    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
853    Qdata = T(Qpoints)
854   
855    qLimits = data['QScaleLim']
856    minQ = np.searchsorted(Qpoints,qLimits[0])
857    maxQ = np.searchsorted(Qpoints,qLimits[1])
858    newdata = []
859    xydata['IofQ'][1][0] = Qpoints
860    xydata['IofQ'][1][1] = Qdata
861    for item in xydata['IofQ'][1]:
862        newdata.append(item[:maxQ])
863    xydata['IofQ'][1] = newdata
864   
865
866    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
867    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
868    Q = xydata['SofQ'][1][0]
869    ruland = Ruland(data['Ruland'],wave,Q,CF)
870#    auxPlot.append([Q,ruland,'Ruland'])     
871    CF *= ruland
872#    auxPlot.append([Q,CF,'CF-Corr'])
873    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
874    xydata['SofQ'][1][1] *= scale
875    xydata['SofQ'][1][1] -= CF
876    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
877    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
878    xydata['SofQ'][1][1] *= scale
879   
880    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
881    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
882    if data['Lorch']:
883        xydata['FofQ'][1][1] *= LorchWeight(Q)
884   
885    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
886    nR = len(xydata['GofR'][1][1])
887    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
888    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
889   
890       
891    return auxPlot
892       
893#testing data
894NeedTestData = True
895def TestData():
896#    global NeedTestData
897    NeedTestData = False
898    global bakType
899    bakType = 'chebyschev'
900    global xdata
901    xdata = np.linspace(4.0,40.0,36000)
902    global parmDict0
903    parmDict0 = {
904        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
905        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
906        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
907        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
908        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
909        'Back0':5.384,'Back1':-0.015,'Back2':.004,
910        }
911    global parmDict1
912    parmDict1 = {
913        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
914        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
915        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
916        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
917        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
918        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
919        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
920        'Back0':36.897,'Back1':-0.508,'Back2':.006,
921        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
922        }
923    global varyList
924    varyList = []
925
926def test0():
927    if NeedTestData: TestData()
928    msg = 'test '
929    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
930    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
931    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
932    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
933    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
934    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
935   
936def test1():
937    if NeedTestData: TestData()
938    time0 = time.time()
939    for i in range(100):
940        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
941    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
942   
943   
944
945if __name__ == '__main__':
946    import GSASIItestplot as plot
947    global plotter
948    plotter = plot.PlotNotebook()
949    test0()
950    print "OK"
951    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.