source: trunk/GSASIIpwd.py @ 1174

Last change on this file since 1174 was 1174, checked in by vondreele, 10 years ago

correct image integration for sample absorption - cylinders only for now
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.9 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2013-12-18 20:55:50 +0000 (Wed, 18 Dec 2013) $
10# $Author: vondreele $
11# $Revision: 1174 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1174 2013-12-18 20:55:50Z vondreele $
14########### SVN repository information ###################
15import sys
16import math
17import time
18
19import numpy as np
20import scipy as sp
21import numpy.linalg as nl
22from numpy.fft import ifft, fft, fftshift
23import scipy.interpolate as si
24import scipy.stats as st
25import scipy.optimize as so
26
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 1174 $")
29import GSASIIlattice as G2lat
30import GSASIIspc as G2spc
31import GSASIIElem as G2elem
32import GSASIIgrid as G2gd
33import GSASIIIO as G2IO
34import GSASIImath as G2mth
35import pypowder as pyd
36
37# trig functions in degrees
38sind = lambda x: math.sin(x*math.pi/180.)
39asind = lambda x: 180.*math.asin(x)/math.pi
40tand = lambda x: math.tan(x*math.pi/180.)
41atand = lambda x: 180.*math.atan(x)/math.pi
42atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
43cosd = lambda x: math.cos(x*math.pi/180.)
44acosd = lambda x: 180.*math.acos(x)/math.pi
45rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
46#numpy versions
47npsind = lambda x: np.sin(x*np.pi/180.)
48npasind = lambda x: 180.*np.arcsin(x)/math.pi
49npcosd = lambda x: np.cos(x*math.pi/180.)
50npacosd = lambda x: 180.*np.arccos(x)/math.pi
51nptand = lambda x: np.tan(x*math.pi/180.)
52npatand = lambda x: 180.*np.arctan(x)/np.pi
53npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
54npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
55npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
56   
57#GSASII pdf calculation routines
58       
59def Transmission(Geometry,Abs,Diam):
60    '''
61    Calculate sample transmission
62
63    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
64    :param float Abs: absorption coeff in cm-1
65    :param float Diam: sample thickness/diameter in mm
66    '''
67    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
68        MuR = Abs*Diam/20.0
69        if MuR <= 3.0:
70            T0 = 16/(3.*math.pi)
71            T1 = -0.045780
72            T2 = -0.02489
73            T3 = 0.003045
74            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
75            if T < -20.:
76                return 2.06e-9
77            else:
78                return math.exp(T)
79        else:
80            T1 = 1.433902
81            T2 = 0.013869+0.337894
82            T3 = 1.933433+1.163198
83            T4 = 0.044365-0.04259
84            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
85            return T/100.
86    elif 'plate' in Geometry:
87        MuR = Abs*Diam/10.
88        return math.exp(-MuR)
89    elif 'Bragg' in Geometry:
90        return 0.0
91
92def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
93    '''Calculate sample absorption
94    :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
95    :param float MuR: absorption coeff * sample thickness/2 or radius
96    :param Tth: 2-theta scattering angle - can be numpy array
97    :param float Phi: flat plate tilt angle - future
98    :param float Psi: flat plate tilt axis - future
99    '''
100   
101    def muRunder3(Sth2):
102        T0 = 16.0/(3.*np.pi)
103        T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
104            0.109561*np.sqrt(Sth2)-26.04556
105        T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
106            1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
107        T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
108        Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
109        return np.exp(Trns)
110   
111    def muRover3(Sth2):
112        T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
113            10.02088*Sth2**3-3.36778*Sth2**4
114        T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
115            (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
116        T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
117            0.13576*np.sqrt(Sth2)+1.163198
118        T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
119        Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
120        return Trns/100.
121       
122    Sth2 = npsind(Tth/2.0)**2
123    Cth2 = 1.-Sth2
124    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
125        return np.where(MuR <3.0,muRunder3(Sth2),muRover3(Sth2))
126    elif 'Bragg' in Geometry:
127        return 1.0
128    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
129        # and only defined for 2theta < 90
130        MuT = 2.*MuR
131        T1 = np.exp(-MuT)
132        T2 = np.exp(-MuT/npcosd(Tth))
133        Tb = MuT-MuT/npcosd(Tth)
134        return (T2-T1)/Tb
135    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
136        MuT = 2.*MuR
137        cth = npcosd(Tth/2.0)
138        return np.exp(-MuT/cth)/cth
139       
140def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
141    'needs a doc string'
142    dA = 0.001
143    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
144    if MuR:
145        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
146        return (AbsP-AbsM)/(2.0*dA)
147    else:
148        return (AbsP-1.)/dA
149       
150def Polarization(Pola,Tth,Azm=0.0):
151    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
152
153    :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
154    :param Azm: azimuthal angle e.g. 0.0 in plane of polarization
155    :param Tth: 2-theta scattering angle - can be numpy array
156      which (if either) of these is "right"?
157    :return: (pola, dpdPola)
158      * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
159        (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
160      * dpdPola: derivative needed for least squares
161
162    """
163    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
164        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
165    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
166    return pola,dpdPola
167   
168def Oblique(ObCoeff,Tth):
169    'currently assumes detector is normal to beam'
170    if ObCoeff:
171        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
172    else:
173        return 1.0
174               
175def Ruland(RulCoff,wave,Q,Compton):
176    'needs a doc string'
177    C = 2.9978e8
178    D = 1.5e-3
179    hmc = 0.024262734687
180    sinth2 = (Q*wave/(4.0*np.pi))**2
181    dlam = (wave**2)*Compton*Q/C
182    dlam_c = 2.0*hmc*sinth2-D*wave**2
183    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
184   
185def LorchWeight(Q):
186    'needs a doc string'
187    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
188           
189def GetAsfMean(ElList,Sthl2):
190    '''Calculate various scattering factor terms for PDF calcs
191
192    :param dict ElList: element dictionary contains scattering factor coefficients, etc.
193    :param np.array Sthl2: numpy array of sin theta/lambda squared values
194    :returns: mean(f^2), mean(f)^2, mean(compton)
195    '''
196    sumNoAtoms = 0.0
197    FF = np.zeros_like(Sthl2)
198    FF2 = np.zeros_like(Sthl2)
199    CF = np.zeros_like(Sthl2)
200    for El in ElList:
201        sumNoAtoms += ElList[El]['FormulaNo']
202    for El in ElList:
203        el = ElList[El]
204        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
205        cf = G2elem.ComptonFac(el,Sthl2)
206        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
207        FF2 += ff2*el['FormulaNo']/sumNoAtoms
208        CF += cf*el['FormulaNo']/sumNoAtoms
209    return FF2,FF**2,CF
210   
211def GetNumDensity(ElList,Vol):
212    'needs a doc string'
213    sumNoAtoms = 0.0
214    for El in ElList:
215        sumNoAtoms += ElList[El]['FormulaNo']
216    return sumNoAtoms/Vol
217           
218def CalcPDF(data,inst,xydata):
219    'needs a doc string'
220    auxPlot = []
221    import copy
222    import scipy.fftpack as ft
223    #subtract backgrounds - if any
224    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
225    if data['Sample Bkg.']['Name']:
226        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
227            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
228    if data['Container']['Name']:
229        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
230        if data['Container Bkg.']['Name']:
231            xycontainer += (xydata['Container Bkg.'][1][1]+
232                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
233        xydata['IofQ'][1][1] += xycontainer
234    #get element data & absorption coeff.
235    ElList = data['ElList']
236    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
237    #Apply angle dependent corrections
238    Tth = xydata['Sample'][1][0]
239    dt = (Tth[1]-Tth[0])
240    MuR = Abs*data['Diam']/20.0
241    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
242    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
243    if data['DetType'] == 'Image plate':
244        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
245    XY = xydata['IofQ'][1]   
246    #convert to Q
247    hc = 12.397639
248    wave = G2mth.getWave(inst)
249    keV = hc/wave
250    minQ = npT2q(Tth[0],wave)
251    maxQ = npT2q(Tth[-1],wave)   
252    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
253    dq = Qpoints[1]-Qpoints[0]
254    XY[0] = npT2q(XY[0],wave)   
255#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
256    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
257    Qdata = T(Qpoints)
258   
259    qLimits = data['QScaleLim']
260    minQ = np.searchsorted(Qpoints,qLimits[0])
261    maxQ = np.searchsorted(Qpoints,qLimits[1])
262    newdata = []
263    xydata['IofQ'][1][0] = Qpoints
264    xydata['IofQ'][1][1] = Qdata
265    for item in xydata['IofQ'][1]:
266        newdata.append(item[:maxQ])
267    xydata['IofQ'][1] = newdata
268   
269
270    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
271    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
272    Q = xydata['SofQ'][1][0]
273    ruland = Ruland(data['Ruland'],wave,Q,CF)
274#    auxPlot.append([Q,ruland,'Ruland'])     
275    CF *= ruland
276#    auxPlot.append([Q,CF,'CF-Corr'])
277    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
278    xydata['SofQ'][1][1] *= scale
279    xydata['SofQ'][1][1] -= CF
280    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
281    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
282    xydata['SofQ'][1][1] *= scale
283   
284    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
285    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
286    if data['Lorch']:
287        xydata['FofQ'][1][1] *= LorchWeight(Q)
288   
289    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
290    nR = len(xydata['GofR'][1][1])
291    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
292    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
293   
294       
295    return auxPlot
296       
297#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
298
299def factorize(num):
300    ''' Provide prime number factors for integer num
301    Returns dictionary of prime factors (keys) & power for each (data)
302    '''
303    factors = {}
304    orig = num
305
306    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
307    i, sqi = 2, 4
308    while sqi <= num:
309        while not num%i:
310            num /= i
311            factors[i] = factors.get(i, 0) + 1
312
313        sqi += 2*i + 1
314        i += 1
315
316    if num != 1 and num != orig:
317        factors[num] = factors.get(num, 0) + 1
318
319    if factors:
320        return factors
321    else:
322        return {num:1}          #a prime number!
323           
324def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
325    ''' Provide list of optimal data sizes for FFT calculations
326
327    :param int nmin: minimum data size >= 1
328    :param int nmax: maximum data size > nmin
329    :param int thresh: maximum prime factor allowed
330    :Returns: list of data sizes where the maximum prime factor is < thresh
331    ''' 
332    plist = []
333    nmin = max(1,nmin)
334    nmax = max(nmin+1,nmax)
335    for p in range(nmin,nmax):
336        if max(factorize(p).keys()) < thresh:
337            plist.append(p)
338    return plist
339
340np.seterr(divide='ignore')
341
342# Normal distribution
343
344# loc = mu, scale = std
345_norm_pdf_C = 1./math.sqrt(2*math.pi)
346class norm_gen(st.rv_continuous):
347    'needs a doc string'
348     
349    def pdf(self,x,*args,**kwds):
350        loc,scale=kwds['loc'],kwds['scale']
351        x = (x-loc)/scale
352        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
353       
354norm = norm_gen(name='norm',longname='A normal',extradoc="""
355
356Normal distribution
357
358The location (loc) keyword specifies the mean.
359The scale (scale) keyword specifies the standard deviation.
360
361normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
362""")
363
364## Cauchy
365
366# median = loc
367
368class cauchy_gen(st.rv_continuous):
369    'needs a doc string'
370
371    def pdf(self,x,*args,**kwds):
372        loc,scale=kwds['loc'],kwds['scale']
373        x = (x-loc)/scale
374        return 1.0/np.pi/(1.0+x*x) / scale
375       
376cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
377
378Cauchy distribution
379
380cauchy.pdf(x) = 1/(pi*(1+x**2))
381
382This is the t distribution with one degree of freedom.
383""")
384   
385   
386#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
387
388
389class fcjde_gen(st.rv_continuous):
390    """
391    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
392    Ref: J. Appl. Cryst. (1994) 27, 892-900.
393
394    :param x: array -1 to 1
395    :param t: 2-theta position of peak
396    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
397      L: sample to detector opening distance
398    :param dx: 2-theta step size in deg
399
400    :returns: for fcj.pdf
401
402     * T = x*dx+t
403     * s = S/L+H/L
404     * if x < 0::
405
406        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
407
408     * if x >= 0: fcj.pdf = 0   
409    """
410    def _pdf(self,x,t,s,dx):
411        T = dx*x+t
412        ax2 = abs(npcosd(T))
413        ax = ax2**2
414        bx = npcosd(t)**2
415        bx = np.where(ax>bx,bx,ax)
416        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
417        fx = np.where(fx > 0.,fx,0.0)
418        return fx
419             
420    def pdf(self,x,*args,**kwds):
421        loc=kwds['loc']
422        return self._pdf(x-loc,*args)
423       
424fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
425               
426def getWidthsCW(pos,sig,gam,shl):
427    'needs a doc string'
428    widths = [np.sqrt(sig)/100.,gam/200.]
429    fwhm = 2.355*widths[0]+2.*widths[1]
430    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
431    fmax = 15.0*fwhm
432    if pos > 90:
433        fmin,fmax = [fmax,fmin]         
434    return widths,fmin,fmax
435   
436def getWidthsTOF(pos,alp,bet,sig,gam):
437    'needs a doc string'
438    lnf = 3.3      # =log(0.001/2)
439    widths = [np.sqrt(sig),gam]
440    fwhm = 2.355*widths[0]+2.*widths[1]
441    fmin = 10.*fwhm*(1.+1./alp)   
442    fmax = 10.*fwhm*(1.+1./bet)
443    return widths,fmin,fmax
444   
445def getFWHM(pos,Inst):
446    'needs a doc string'
447    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
448    sigTOF = lambda dsp,S0,S1,Sq:  S0+S1*dsp**2+Sq*dsp
449    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
450    gamTOF = lambda dsp,X,Y: X*dsp+Y*dsp**2
451    if 'C' in Inst['Type'][0]:
452        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
453        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])*100.
454    else:
455        dsp = pos/Inst['difC'][0]
456        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-q'][1])
457        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
458    return getgamFW(g,s)
459   
460def getgamFW(g,s):
461    'needs a doc string'
462    gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
463    return gamFW(g,s)
464               
465def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
466    'needs a doc string'
467    DX = xdata[1]-xdata[0]
468    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
469    x = np.linspace(pos-fmin,pos+fmin,256)
470    dx = x[1]-x[0]
471    Norm = norm.pdf(x,loc=pos,scale=widths[0])
472    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
473    arg = [pos,shl/57.2958,dx,]
474    FCJ = fcjde.pdf(x,*arg,loc=pos)
475    if len(np.nonzero(FCJ)[0])>5:
476        z = np.column_stack([Norm,Cauchy,FCJ]).T
477        Z = fft(z)
478        Df = ifft(Z.prod(axis=0)).real
479    else:
480        z = np.column_stack([Norm,Cauchy]).T
481        Z = fft(z)
482        Df = fftshift(ifft(Z.prod(axis=0))).real
483    Df /= np.sum(Df)
484    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
485    return intens*Df(xdata)*DX/dx
486
487def getBackground(pfx,parmDict,bakType,xdata):
488    'needs a doc string'
489    yb = np.zeros_like(xdata)
490    nBak = 0
491    cw = np.diff(xdata)
492    cw = np.append(cw,cw[-1])
493    while True:
494        key = pfx+'Back:'+str(nBak)
495        if key in parmDict:
496            nBak += 1
497        else:
498            break
499    if bakType in ['chebyschev','cosine']:
500        dt = xdata[-1]-xdata[0]   
501        for iBak in range(nBak):
502            key = pfx+'Back:'+str(iBak)
503            if bakType == 'chebyschev':
504                yb += parmDict[key]*(2.*(xdata-xdata[0])/dt-1.)**iBak
505            elif bakType == 'cosine':
506                yb += parmDict[key]*npcosd(xdata*iBak)
507    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
508        if nBak == 1:
509            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
510        elif nBak == 2:
511            dX = xdata[-1]-xdata[0]
512            T2 = (xdata-xdata[0])/dX
513            T1 = 1.0-T2
514            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
515        else:
516            if bakType == 'lin interpolate':
517                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
518            elif bakType == 'inv interpolate':
519                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
520            elif bakType == 'log interpolate':
521                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
522            bakPos[0] = xdata[0]
523            bakPos[-1] = xdata[-1]
524            bakVals = np.zeros(nBak)
525            for i in range(nBak):
526                bakVals[i] = parmDict[pfx+'Back:'+str(i)]
527            bakInt = si.interp1d(bakPos,bakVals,'linear')
528            yb = bakInt(xdata)
529    if 'difC' in parmDict:
530        ff = 1.
531    else:       
532        try:
533            wave = parmDict[pfx+'Lam']
534        except KeyError:
535            wave = parmDict[pfx+'Lam1']
536        q = 4.0*np.pi*npsind(xdata/2.0)/wave
537        SQ = (q/(4.*np.pi))**2
538        FF = G2elem.GetFormFactorCoeff('Si')[0]
539        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
540    iD = 0       
541    while True:
542        try:
543            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
544            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
545            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
546            yb += ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
547            iD += 1       
548        except KeyError:
549            break
550    iD = 0
551    while True:
552        try:
553            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
554            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
555            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
556            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
557            shl = 0.002
558            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
559            iBeg = np.searchsorted(xdata,pkP-fmin)
560            iFin = np.searchsorted(xdata,pkP+fmax)
561            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
562            iD += 1       
563        except KeyError:
564            break
565        except ValueError:
566            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
567            break       
568    return yb
569   
570def getBackgroundDerv(pfx,parmDict,bakType,xdata):
571    'needs a doc string'
572    nBak = 0
573    while True:
574        key = pfx+'Back:'+str(nBak)
575        if key in parmDict:
576            nBak += 1
577        else:
578            break
579    dydb = np.zeros(shape=(nBak,len(xdata)))
580    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
581    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
582    cw = np.diff(xdata)
583    cw = np.append(cw,cw[-1])
584
585    if bakType in ['chebyschev','cosine']:
586        dt = xdata[-1]-xdata[0]   
587        for iBak in range(nBak):   
588            if bakType == 'chebyschev':
589                dydb[iBak] = (2.*(xdata-xdata[0])/dt-1.)**iBak
590            elif bakType == 'cosine':
591                dydb[iBak] = npcosd(xdata*iBak)
592    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
593        if nBak == 1:
594            dydb[0] = np.ones_like(xdata)
595        elif nBak == 2:
596            dX = xdata[-1]-xdata[0]
597            T2 = (xdata-xdata[0])/dX
598            T1 = 1.0-T2
599            dydb = [T1,T2]
600        else:
601            if bakType == 'lin interpolate':
602                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
603            elif bakType == 'inv interpolate':
604                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
605            elif bakType == 'log interpolate':
606                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
607            bakPos[0] = xdata[0]
608            bakPos[-1] = xdata[-1]
609            for i,pos in enumerate(bakPos):
610                if i == 0:
611                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
612                elif i == len(bakPos)-1:
613                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
614                else:
615                    dydb[i] = np.where(xdata>bakPos[i],
616                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
617                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
618    if 'difC' in parmDict:
619        ff = 1.
620    else:
621        try:
622            wave = parmDict[pfx+'Lam']
623        except KeyError:
624            wave = parmDict[pfx+'Lam1']
625        q = 4.0*np.pi*npsind(xdata/2.0)/wave
626        SQ = (q/(4*np.pi))**2
627        FF = G2elem.GetFormFactorCoeff('Si')[0]
628        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
629    iD = 0       
630    while True:
631        try:
632            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
633            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
634            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
635            sqr = np.sin(q*dbR)/(q*dbR)
636            cqr = np.cos(q*dbR)
637            temp = np.exp(-dbU*q**2)
638            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
639            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
640            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
641            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
642        except KeyError:
643            break
644    iD = 0
645    while True:
646        try:
647            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
648            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
649            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
650            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
651            shl = 0.002
652            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
653            iBeg = np.searchsorted(xdata,pkP-fmin)
654            iFin = np.searchsorted(xdata,pkP+fmax)
655            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
656            dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
657            dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
658            dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
659            dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
660            iD += 1       
661        except KeyError:
662            break
663        except ValueError:
664            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
665            break       
666    return dydb,dyddb,dydpk
667
668#use old fortran routine
669def getFCJVoigt3(pos,sig,gam,shl,xdata):
670    'needs a doc string'
671   
672    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
673#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
674    Df /= np.sum(Df)
675    return Df
676
677def getdFCJVoigt3(pos,sig,gam,shl,xdata):
678    'needs a doc string'
679   
680    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
681#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
682    sumDf = np.sum(Df)
683    return Df,dFdp,dFds,dFdg,dFdsh
684
685def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
686    'needs a doc string'
687    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
688    Df /= np.sum(Df)
689    return Df 
690   
691def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
692    'needs a doc string'
693    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
694    sumDf = np.sum(Df)
695    return Df,dFdp,dFda,dFdb,dFds,dFdg   
696
697def ellipseSize(H,Sij,GB):
698    'needs a doc string'
699    HX = np.inner(H.T,GB)
700    lenHX = np.sqrt(np.sum(HX**2))
701    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
702    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
703    lenR = np.sqrt(np.sum(R**2))
704    return lenR
705
706def ellipseSizeDerv(H,Sij,GB):
707    'needs a doc string'
708    lenR = ellipseSize(H,Sij,GB)
709    delt = 0.001
710    dRdS = np.zeros(6)
711    for i in range(6):
712        dSij = Sij[:]
713        dSij[i] += delt
714        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
715    return lenR,dRdS
716
717def getHKLpeak(dmin,SGData,A):
718    'needs a doc string'
719    HKL = G2lat.GenHLaue(dmin,SGData,A)       
720    HKLs = []
721    for h,k,l,d in HKL:
722        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
723        if not ext:
724            HKLs.append([h,k,l,d,-1])
725    return HKLs
726
727def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
728    'needs a doc string'
729   
730    yb = getBackground('',parmDict,bakType,xdata)
731    yc = np.zeros_like(yb)
732    cw = np.diff(xdata)
733    cw = np.append(cw,cw[-1])
734    if 'C' in dataType:
735        shl = max(parmDict['SH/L'],0.002)
736        Ka2 = False
737        if 'Lam1' in parmDict.keys():
738            Ka2 = True
739            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
740            kRatio = parmDict['I(L2)/I(L1)']
741        iPeak = 0
742        while True:
743            try:
744                pos = parmDict['pos'+str(iPeak)]
745                theta = (pos-parmDict['Zero'])/2.0
746                intens = parmDict['int'+str(iPeak)]
747                sigName = 'sig'+str(iPeak)
748                if sigName in varyList:
749                    sig = parmDict[sigName]
750                else:
751                    sig = G2mth.getCWsig(parmDict,theta)
752                sig = max(sig,0.001)          #avoid neg sigma
753                gamName = 'gam'+str(iPeak)
754                if gamName in varyList:
755                    gam = parmDict[gamName]
756                else:
757                    gam = G2mth.getCWgam(parmDict,theta)
758                gam = max(gam,0.001)             #avoid neg gamma
759                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
760                iBeg = np.searchsorted(xdata,pos-fmin)
761                iFin = np.searchsorted(xdata,pos+fmin)
762                if not iBeg+iFin:       #peak below low limit
763                    iPeak += 1
764                    continue
765                elif not iBeg-iFin:     #peak above high limit
766                    return yb+yc
767                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
768                if Ka2:
769                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
770                    iBeg = np.searchsorted(xdata,pos2-fmin)
771                    iFin = np.searchsorted(xdata,pos2+fmin)
772                    if iBeg-iFin:
773                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
774                iPeak += 1
775            except KeyError:        #no more peaks to process
776                return yb+yc
777    else:
778        Pdabc = parmDict['Pdabc']
779        difC = parmDict['difC']
780        iPeak = 0
781        while True:
782            try:
783                pos = parmDict['pos'+str(iPeak)]               
784                tof = pos-parmDict['Zero']
785                dsp = tof/difC
786                intens = parmDict['int'+str(iPeak)]
787                alpName = 'alp'+str(iPeak)
788                if alpName in varyList:
789                    alp = parmDict[alpName]
790                else:
791                    if len(Pdabc):
792                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
793                    else:
794                        alp = G2mth.getTOFalpha(parmDict,dsp)
795                betName = 'bet'+str(iPeak)
796                if betName in varyList:
797                    bet = parmDict[betName]
798                else:
799                    if len(Pdabc):
800                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
801                    else:
802                        bet = G2mth.getTOFbeta(parmDict,dsp)
803                sigName = 'sig'+str(iPeak)
804                if sigName in varyList:
805                    sig = parmDict[sigName]
806                else:
807                    sig = G2mth.getTOFsig(parmDict,dsp)
808                gamName = 'gam'+str(iPeak)
809                if gamName in varyList:
810                    gam = parmDict[gamName]
811                else:
812                    gam = G2mth.getTOFgamma(parmDict,dsp)
813                gam = max(gam,0.001)             #avoid neg gamma
814                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
815                iBeg = np.searchsorted(xdata,pos-fmin)
816                iFin = np.searchsorted(xdata,pos+fmax)
817                lenX = len(xdata)
818                if not iBeg:
819                    iFin = np.searchsorted(xdata,pos+fmax)
820                elif iBeg == lenX:
821                    iFin = iBeg
822                else:
823                    iFin = np.searchsorted(xdata,pos+fmax)
824                if not iBeg+iFin:       #peak below low limit
825                    iPeak += 1
826                    continue
827                elif not iBeg-iFin:     #peak above high limit
828                    return yb+yc
829                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
830                iPeak += 1
831            except KeyError:        #no more peaks to process
832                return yb+yc
833           
834def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
835    'needs a doc string'
836# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
837    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
838    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
839    if 'Back:0' in varyList:            #background derivs are in front if present
840        dMdv[0:len(dMdb)] = dMdb
841    names = ['DebyeA','DebyeR','DebyeU']
842    for name in varyList:
843        if 'Debye' in name:
844            parm,id = name.split(':')
845            ip = names.index(parm)
846            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
847    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
848    for name in varyList:
849        if 'BkPk' in name:
850            parm,id = name.split(':')
851            ip = names.index(parm)
852            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
853    cw = np.diff(xdata)
854    cw = np.append(cw,cw[-1])
855    if 'C' in dataType:
856        shl = max(parmDict['SH/L'],0.002)
857        Ka2 = False
858        if 'Lam1' in parmDict.keys():
859            Ka2 = True
860            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
861            kRatio = parmDict['I(L2)/I(L1)']
862        iPeak = 0
863        while True:
864            try:
865                pos = parmDict['pos'+str(iPeak)]
866                theta = (pos-parmDict['Zero'])/2.0
867                intens = parmDict['int'+str(iPeak)]
868                sigName = 'sig'+str(iPeak)
869                tanth = tand(theta)
870                costh = cosd(theta)
871                if sigName in varyList:
872                    sig = parmDict[sigName]
873                else:
874                    sig = G2mth.getCWsig(parmDict,theta)
875                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(theta)
876                sig = max(sig,0.001)          #avoid neg sigma
877                gamName = 'gam'+str(iPeak)
878                if gamName in varyList:
879                    gam = parmDict[gamName]
880                else:
881                    gam = G2mth.getCWgam(parmDict,theta)
882                    dgdX,dgdY = G2mth.getCWgamDeriv(theta)
883                gam = max(gam,0.001)             #avoid neg gamma
884                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
885                iBeg = np.searchsorted(xdata,pos-fmin)
886                iFin = np.searchsorted(xdata,pos+fmin)
887                if not iBeg+iFin:       #peak below low limit
888                    iPeak += 1
889                    continue
890                elif not iBeg-iFin:     #peak above high limit
891                    break
892                dMdpk = np.zeros(shape=(6,len(xdata)))
893                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
894                for i in range(1,5):
895                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
896                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
897                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
898                if Ka2:
899                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
900                    iBeg = np.searchsorted(xdata,pos2-fmin)
901                    iFin = np.searchsorted(xdata,pos2+fmin)
902                    if iBeg-iFin:
903                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
904                        for i in range(1,5):
905                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
906                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
907                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
908                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
909                for parmName in ['pos','int','sig','gam']:
910                    try:
911                        idx = varyList.index(parmName+str(iPeak))
912                        dMdv[idx] = dervDict[parmName]
913                    except ValueError:
914                        pass
915                if 'U' in varyList:
916                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
917                if 'V' in varyList:
918                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
919                if 'W' in varyList:
920                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
921                if 'X' in varyList:
922                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
923                if 'Y' in varyList:
924                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
925                if 'SH/L' in varyList:
926                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
927                if 'I(L2)/I(L1)' in varyList:
928                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
929                iPeak += 1
930            except KeyError:        #no more peaks to process
931                break
932    else:
933        Pdabc = parmDict['Pdabc']
934        difC = parmDict['difC']
935        iPeak = 0
936        while True:
937            try:
938                pos = parmDict['pos'+str(iPeak)]               
939                tof = pos-parmDict['Zero']
940                dsp = tof/difC
941                intens = parmDict['int'+str(iPeak)]
942                alpName = 'alp'+str(iPeak)
943                if alpName in varyList:
944                    alp = parmDict[alpName]
945                else:
946                    if len(Pdabc):
947                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
948                    else:
949                        alp = G2mth.getTOFalpha(parmDict,dsp)
950                        dada0 = G2mth.getTOFalphaDeriv(dsp)
951                betName = 'bet'+str(iPeak)
952                if betName in varyList:
953                    bet = parmDict[betName]
954                else:
955                    if len(Pdabc):
956                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
957                    else:
958                        bet = G2mth.getTOFbeta(parmDict,dsp)
959                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
960                sigName = 'sig'+str(iPeak)
961                if sigName in varyList:
962                    sig = parmDict[sigName]
963                else:
964                    sig = G2mth.getTOFsig(parmDict,dsp)
965                    dsds0,dsds1,dsds2 = G2mth.getTOFsigDeriv(dsp)
966                gamName = 'gam'+str(iPeak)
967                if gamName in varyList:
968                    gam = parmDict[gamName]
969                else:
970                    gam = G2mth.getTOFgamma(parmDict,dsp)
971                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
972                gam = max(gam,0.001)             #avoid neg gamma
973                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
974                iBeg = np.searchsorted(xdata,pos-fmin)
975                lenX = len(xdata)
976                if not iBeg:
977                    iFin = np.searchsorted(xdata,pos+fmax)
978                elif iBeg == lenX:
979                    iFin = iBeg
980                else:
981                    iFin = np.searchsorted(xdata,pos+fmax)
982                if not iBeg+iFin:       #peak below low limit
983                    iPeak += 1
984                    continue
985                elif not iBeg-iFin:     #peak above high limit
986                    break
987                dMdpk = np.zeros(shape=(7,len(xdata)))
988                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
989                for i in range(1,6):
990                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
991                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
992                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
993                for parmName in ['pos','int','alp','bet','sig','gam']:
994                    try:
995                        idx = varyList.index(parmName+str(iPeak))
996                        dMdv[idx] = dervDict[parmName]
997                    except ValueError:
998                        pass
999                if 'alpha' in varyList:
1000                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1001                if 'beta-0' in varyList:
1002                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1003                if 'beta-1' in varyList:
1004                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1005                if 'beta-q' in varyList:
1006                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1007                if 'sig-0' in varyList:
1008                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1009                if 'sig-1' in varyList:
1010                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1011                if 'sig-q' in varyList:
1012                    dMdv[varyList.index('sig-q')] += dsds2*dervDict['sig']
1013                if 'X' in varyList:
1014                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1015                if 'Y' in varyList:
1016                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1017                iPeak += 1
1018            except KeyError:        #no more peaks to process
1019                break
1020    return dMdv
1021       
1022def Dict2Values(parmdict, varylist):
1023    '''Use before call to leastsq to setup list of values for the parameters
1024    in parmdict, as selected by key in varylist'''
1025    return [parmdict[key] for key in varylist] 
1026   
1027def Values2Dict(parmdict, varylist, values):
1028    ''' Use after call to leastsq to update the parameter dictionary with
1029    values corresponding to keys in varylist'''
1030    parmdict.update(zip(varylist,values))
1031   
1032def SetBackgroundParms(Background):
1033    'needs a doc string'
1034    if len(Background) == 1:            # fix up old backgrounds
1035        BackGround.append({'nDebye':0,'debyeTerms':[]})
1036    bakType,bakFlag = Background[0][:2]
1037    backVals = Background[0][3:]
1038    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1039    Debye = Background[1]           #also has background peaks stuff
1040    backDict = dict(zip(backNames,backVals))
1041    backVary = []
1042    if bakFlag:
1043        backVary = backNames
1044
1045    backDict['nDebye'] = Debye['nDebye']
1046    debyeDict = {}
1047    debyeList = []
1048    for i in range(Debye['nDebye']):
1049        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1050        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1051        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1052    debyeVary = []
1053    for item in debyeList:
1054        if item[1]:
1055            debyeVary.append(item[0])
1056    backDict.update(debyeDict)
1057    backVary += debyeVary
1058
1059    backDict['nPeaks'] = Debye['nPeaks']
1060    peaksDict = {}
1061    peaksList = []
1062    for i in range(Debye['nPeaks']):
1063        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1064        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1065        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1066    peaksVary = []
1067    for item in peaksList:
1068        if item[1]:
1069            peaksVary.append(item[0])
1070    backDict.update(peaksDict)
1071    backVary += peaksVary   
1072    return bakType,backDict,backVary
1073           
1074def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1075    'needs a doc string'
1076   
1077       
1078    def GetBackgroundParms(parmList,Background):
1079        iBak = 0
1080        while True:
1081            try:
1082                bakName = 'Back:'+str(iBak)
1083                Background[0][iBak+3] = parmList[bakName]
1084                iBak += 1
1085            except KeyError:
1086                break
1087        iDb = 0
1088        while True:
1089            names = ['DebyeA:','DebyeR:','DebyeU:']
1090            try:
1091                for i,name in enumerate(names):
1092                    val = parmList[name+str(iDb)]
1093                    Background[1]['debyeTerms'][iDb][2*i] = val
1094                iDb += 1
1095            except KeyError:
1096                break
1097        iDb = 0
1098        while True:
1099            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1100            try:
1101                for i,name in enumerate(names):
1102                    val = parmList[name+str(iDb)]
1103                    Background[1]['peaksList'][iDb][2*i] = val
1104                iDb += 1
1105            except KeyError:
1106                break
1107               
1108    def BackgroundPrint(Background,sigDict):
1109        if Background[0][1]:
1110            print 'Background coefficients for',Background[0][0],'function'
1111            ptfmt = "%12.5f"
1112            ptstr =  'values:'
1113            sigstr = 'esds  :'
1114            for i,back in enumerate(Background[0][3:]):
1115                ptstr += ptfmt % (back)
1116                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1117            print ptstr
1118            print sigstr
1119        else:
1120            print 'Background not refined'
1121        if Background[1]['nDebye']:
1122            parms = ['DebyeA','DebyeR','DebyeU']
1123            print 'Debye diffuse scattering coefficients'
1124            ptfmt = "%12.5f"
1125            names =   'names :'
1126            ptstr =  'values:'
1127            sigstr = 'esds  :'
1128            for item in sigDict:
1129                if 'Debye' in item:
1130                    names += '%12s'%(item)
1131                    sigstr += ptfmt%(sigDict[item])
1132                    parm,id = item.split(':')
1133                    ip = parms.index(parm)
1134                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1135            print names
1136            print ptstr
1137            print sigstr
1138        if Background[1]['nPeaks']:
1139            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1140            print 'Peaks in background coefficients'
1141            ptfmt = "%12.5f"
1142            names =   'names :'
1143            ptstr =  'values:'
1144            sigstr = 'esds  :'
1145            for item in sigDict:
1146                if 'BkPk' in item:
1147                    names += '%12s'%(item)
1148                    sigstr += ptfmt%(sigDict[item])
1149                    parm,id = item.split(':')
1150                    ip = parms.index(parm)
1151                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1152            print names
1153            print ptstr
1154            print sigstr
1155                           
1156    def SetInstParms(Inst):
1157        dataType = Inst['Type'][0]
1158        insVary = []
1159        insNames = []
1160        insVals = []
1161        for parm in Inst:
1162            insNames.append(parm)
1163            insVals.append(Inst[parm][1])
1164            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1165                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',] and Inst[parm][2]:
1166                    insVary.append(parm)
1167        instDict = dict(zip(insNames,insVals))
1168        instDict['X'] = max(instDict['X'],0.01)
1169        instDict['Y'] = max(instDict['Y'],0.01)
1170        if 'SH/L' in instDict:
1171            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1172        return dataType,instDict,insVary
1173       
1174    def GetInstParms(parmDict,Inst,varyList,Peaks):
1175        for name in Inst:
1176            Inst[name][1] = parmDict[name]
1177        iPeak = 0
1178        while True:
1179            try:
1180                sigName = 'sig'+str(iPeak)
1181                pos = parmDict['pos'+str(iPeak)]
1182                if sigName not in varyList:
1183                    if 'C' in Inst['Type'][0]:
1184                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1185                    else:
1186                        dsp = pos/Inst['difC'][1]
1187                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1188                gamName = 'gam'+str(iPeak)
1189                if gamName not in varyList:
1190                    if 'C' in Inst['Type'][0]:
1191                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1192                    else:
1193                        dsp = pos/Inst['difC'][1]
1194                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1195                iPeak += 1
1196            except KeyError:
1197                break
1198       
1199    def InstPrint(Inst,sigDict):
1200        print 'Instrument Parameters:'
1201        ptfmt = "%12.6f"
1202        ptlbls = 'names :'
1203        ptstr =  'values:'
1204        sigstr = 'esds  :'
1205        for parm in Inst:
1206            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1207                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',]:
1208                ptlbls += "%s" % (parm.center(12))
1209                ptstr += ptfmt % (Inst[parm][1])
1210                if parm in sigDict:
1211                    sigstr += ptfmt % (sigDict[parm])
1212                else:
1213                    sigstr += 12*' '
1214        print ptlbls
1215        print ptstr
1216        print sigstr
1217
1218    def SetPeaksParms(dataType,Peaks):
1219        peakNames = []
1220        peakVary = []
1221        peakVals = []
1222        if 'C' in dataType:
1223            names = ['pos','int','sig','gam']
1224        else:
1225            names = ['pos','int','alp','bet','sig','gam']
1226        for i,peak in enumerate(Peaks):
1227            for j,name in enumerate(names):
1228                peakVals.append(peak[2*j])
1229                parName = name+str(i)
1230                peakNames.append(parName)
1231                if peak[2*j+1]:
1232                    peakVary.append(parName)
1233        return dict(zip(peakNames,peakVals)),peakVary
1234               
1235    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1236        if 'C' in Inst['Type'][0]:
1237            names = ['pos','int','sig','gam']
1238        else:
1239            names = ['pos','int','alp','bet','sig','gam']
1240        for i,peak in enumerate(Peaks):
1241            pos = parmDict['pos'+str(i)]
1242            if 'difC' in Inst:
1243                dsp = pos/Inst['difC'][1]
1244            for j in range(len(names)):
1245                parName = names[j]+str(i)
1246                if parName in varyList:
1247                    peak[2*j] = parmDict[parName]
1248                elif 'alpha' in parName:
1249                    peak[2*j] = parmDict['alpha']/dsp
1250                elif 'beta' in parName:
1251                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1252                elif 'sig' in parName:
1253                    if 'C' in Inst['Type'][0]:
1254                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1255                    else:
1256                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1257                elif 'gam' in parName:
1258                    if 'C' in Inst['Type'][0]:
1259                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1260                    else:
1261                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1262                       
1263    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1264        print 'Peak coefficients:'
1265        if 'C' in dataType:
1266            names = ['pos','int','sig','gam']
1267        else:
1268            names = ['pos','int','alp','bet','sig','gam']           
1269        head = 13*' '
1270        for name in names:
1271            head += name.center(10)+'esd'.center(10)
1272        print head
1273        if 'C' in dataType:
1274            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1275        else:
1276            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1277        for i,peak in enumerate(Peaks):
1278            ptstr =  ':'
1279            for j in range(len(names)):
1280                name = names[j]
1281                parName = name+str(i)
1282                ptstr += ptfmt[name] % (parmDict[parName])
1283                if parName in varyList:
1284#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1285                    ptstr += ptfmt[name] % (sigDict[parName])
1286                else:
1287#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1288                    ptstr += 10*' '
1289            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1290               
1291    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1292        parmdict.update(zip(varylist,values))
1293        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1294           
1295    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1296        parmdict.update(zip(varylist,values))
1297        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1298        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1299        if dlg:
1300            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1301            if not GoOn:
1302                return -M           #abort!!
1303        return M
1304       
1305    if controls:
1306        Ftol = controls['min dM/M']
1307        derivType = controls['deriv type']
1308    else:
1309        Ftol = 0.0001
1310        derivType = 'analytic'
1311    if oneCycle:
1312        Ftol = 1.0
1313    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1314    yc *= 0.                            #set calcd ones to zero
1315    yb *= 0.
1316    yd *= 0.
1317    xBeg = np.searchsorted(x,Limits[0])
1318    xFin = np.searchsorted(x,Limits[1])
1319    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1320    dataType,insDict,insVary = SetInstParms(Inst)
1321    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1322    parmDict = {}
1323    parmDict.update(bakDict)
1324    parmDict.update(insDict)
1325    parmDict.update(peakDict)
1326    parmDict['Pdabc'] = []      #dummy Pdabc
1327    parmDict.update(Inst2)      #put in real one if there
1328    varyList = bakVary+insVary+peakVary
1329    while True:
1330        begin = time.time()
1331        values =  np.array(Dict2Values(parmDict, varyList))
1332        if FitPgm == 'LSQ':
1333            try:
1334                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1335                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1336                ncyc = int(result[2]['nfev']/2)
1337            finally:
1338                dlg.Destroy()
1339            runtime = time.time()-begin   
1340            chisq = np.sum(result[2]['fvec']**2)
1341            Values2Dict(parmDict, varyList, result[0])
1342            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1343            GOF = chisq/(xFin-xBeg-len(varyList))
1344            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1345            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1346            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1347            try:
1348                sig = np.sqrt(np.diag(result[1])*GOF)
1349                if np.any(np.isnan(sig)):
1350                    print '*** Least squares aborted - some invalid esds possible ***'
1351                break                   #refinement succeeded - finish up!
1352            except ValueError:          #result[1] is None on singular matrix
1353                print '**** Refinement failed - singular matrix ****'
1354                Ipvt = result[2]['ipvt']
1355                for i,ipvt in enumerate(Ipvt):
1356                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1357                        print 'Removing parameter: ',varyList[ipvt-1]
1358                        del(varyList[ipvt-1])
1359                        break
1360        elif FitPgm == 'BFGS':
1361            print 'Other program here'
1362            return
1363       
1364    sigDict = dict(zip(varyList,sig))
1365    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1366    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1367    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1368    GetBackgroundParms(parmDict,Background)
1369    BackgroundPrint(Background,sigDict)
1370    GetInstParms(parmDict,Inst,varyList,Peaks)
1371    InstPrint(Inst,sigDict)
1372    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1373    PeaksPrint(dataType,parmDict,sigDict,varyList)
1374
1375def calcIncident(Iparm,xdata):
1376    'needs a doc string'
1377
1378    def IfunAdv(Iparm,xdata):
1379        Itype = Iparm['Itype']
1380        Icoef = Iparm['Icoeff']
1381        DYI = np.ones((12,xdata.shape[0]))
1382        YI = np.ones_like(xdata)*Icoef[0]
1383       
1384        x = xdata/1000.                 #expressions are in ms
1385        if Itype == 'Exponential':
1386            for i in range(1,10,2):
1387                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1388                YI += Icoef[i]*Eterm
1389                DYI[i] *= Eterm
1390                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1391        elif 'Maxwell'in Itype:
1392            Eterm = np.exp(-Icoef[2]/x**2)
1393            DYI[1] = Eterm/x**5
1394            DYI[2] = -Icoef[1]*DYI[1]/x**2
1395            YI += (Icoef[1]*Eterm/x**5)
1396            if 'Exponential' in Itype:
1397                for i in range(3,12,2):
1398                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1399                    YI += Icoef[i]*Eterm
1400                    DYI[i] *= Eterm
1401                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1402            else:   #Chebyschev
1403                T = (2./x)-1.
1404                Ccof = np.ones((12,xdata.shape[0]))
1405                Ccof[1] = T
1406                for i in range(2,12):
1407                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1408                for i in range(1,10):
1409                    YI += Ccof[i]*Icoef[i+2]
1410                    DYI[i+2] =Ccof[i]
1411        return YI,DYI
1412       
1413    Iesd = np.array(Iparm['Iesd'])
1414    Icovar = Iparm['Icovar']
1415    YI,DYI = IfunAdv(Iparm,xdata)
1416    YI = np.where(YI>0,YI,1.)
1417    WYI = np.zeros_like(xdata)
1418    vcov = np.zeros((12,12))
1419    k = 0
1420    for i in range(12):
1421        for j in range(i,12):
1422            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1423            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1424            k += 1
1425    M = np.inner(vcov,DYI.T)
1426    WYI = np.sum(M*DYI,axis=0)
1427    WYI = np.where(WYI>0.,WYI,0.)
1428    return YI,WYI
1429   
1430#testing data
1431NeedTestData = True
1432def TestData():
1433    'needs a doc string'
1434#    global NeedTestData
1435    NeedTestData = False
1436    global bakType
1437    bakType = 'chebyschev'
1438    global xdata
1439    xdata = np.linspace(4.0,40.0,36000)
1440    global parmDict0
1441    parmDict0 = {
1442        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1443        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1444        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1445        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1446        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1447        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1448        }
1449    global parmDict1
1450    parmDict1 = {
1451        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1452        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1453        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1454        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1455        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1456        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1457        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1458        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1459        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1460        }
1461    global parmDict2
1462    parmDict2 = {
1463        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1464        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1465        'Back0':5.,'Back1':-0.02,'Back2':.004,
1466#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1467        }
1468    global varyList
1469    varyList = []
1470
1471def test0():
1472    if NeedTestData: TestData()
1473    msg = 'test '
1474    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1475    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1476    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1477    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1478    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1479    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1480   
1481def test1():
1482    if NeedTestData: TestData()
1483    time0 = time.time()
1484    for i in range(100):
1485        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1486    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1487   
1488def test2(name,delt):
1489    if NeedTestData: TestData()
1490    varyList = [name,]
1491    xdata = np.linspace(5.6,5.8,400)
1492    hplot = plotter.add('derivatives test for '+name).gca()
1493    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1494    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1495    parmDict2[name] += delt
1496    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1497    hplot.plot(xdata,(y1-y0)/delt,'r+')
1498   
1499def test3(name,delt):
1500    if NeedTestData: TestData()
1501    names = ['pos','sig','gam','shl']
1502    idx = names.index(name)
1503    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1504    xdata = np.linspace(5.6,5.8,800)
1505    dx = xdata[1]-xdata[0]
1506    hplot = plotter.add('derivatives test for '+name).gca()
1507    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1508    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1509    myDict[name] += delt
1510    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1511    hplot.plot(xdata,(y1-y0)/delt,'r+')
1512
1513if __name__ == '__main__':
1514    import GSASIItestplot as plot
1515    global plotter
1516    plotter = plot.PlotNotebook()
1517#    test0()
1518#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1519    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1520        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1521        test2(name,shft)
1522    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1523        test3(name,shft)
1524    print "OK"
1525    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.