source: trunk/GSASIIpwd.py @ 1009

Last change on this file since 1009 was 1009, checked in by vondreele, 9 years ago

auto peak search for TOF - no duplicates either
small amendment to peak fit tutorial

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