source: trunk/GSASIIpwd.py @ 939

Last change on this file since 939 was 939, checked in by toby, 9 years ago

fix & cleanup unit tests; add/change doc strings for sphinx; add all G2 py files to sphinx

  • 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-06-03 04:07:35 +0000 (Mon, 03 Jun 2013) $
10# $Author: toby $
11# $Revision: 939 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 939 2013-06-03 04:07:35Z toby $
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: 939 $")
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(TTh,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    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
445    gamFW = lambda s,g: math.exp(math.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.)
446    s = sig(TTh/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
447    g = gam(TTh/2.,Inst['X'][1],Inst['Y'][1])*100.
448    return gamFW(g,s)   
449               
450def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
451    'needs a doc string'
452    DX = xdata[1]-xdata[0]
453    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
454    x = np.linspace(pos-fmin,pos+fmin,256)
455    dx = x[1]-x[0]
456    Norm = norm.pdf(x,loc=pos,scale=widths[0])
457    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
458    arg = [pos,shl/57.2958,dx,]
459    FCJ = fcjde.pdf(x,*arg,loc=pos)
460    if len(np.nonzero(FCJ)[0])>5:
461        z = np.column_stack([Norm,Cauchy,FCJ]).T
462        Z = fft(z)
463        Df = ifft(Z.prod(axis=0)).real
464    else:
465        z = np.column_stack([Norm,Cauchy]).T
466        Z = fft(z)
467        Df = fftshift(ifft(Z.prod(axis=0))).real
468    Df /= np.sum(Df)
469    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
470    return intens*Df(xdata)*DX/dx
471
472def getBackground(pfx,parmDict,bakType,xdata):
473    'needs a doc string'
474    yb = np.zeros_like(xdata)
475    nBak = 0
476    cw = np.diff(xdata)
477    cw = np.append(cw,cw[-1])
478    while True:
479        key = pfx+'Back:'+str(nBak)
480        if key in parmDict:
481            nBak += 1
482        else:
483            break
484    if bakType in ['chebyschev','cosine']:
485        for iBak in range(nBak):   
486            key = pfx+'Back:'+str(iBak)
487            if bakType == 'chebyschev':
488                yb += parmDict[key]*(xdata-xdata[0])**iBak
489            elif bakType == 'cosine':
490                yb += parmDict[key]*npcosd(xdata*iBak)
491    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
492        if nBak == 1:
493            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
494        elif nBak == 2:
495            dX = xdata[-1]-xdata[0]
496            T2 = (xdata-xdata[0])/dX
497            T1 = 1.0-T2
498            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
499        else:
500            if bakType == 'lin interpolate':
501                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
502            elif bakType == 'inv interpolate':
503                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
504            elif bakType == 'log interpolate':
505                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
506            bakPos[0] = xdata[0]
507            bakPos[-1] = xdata[-1]
508            bakVals = np.zeros(nBak)
509            for i in range(nBak):
510                bakVals[i] = parmDict[pfx+'Back:'+str(i)]
511            bakInt = si.interp1d(bakPos,bakVals,'linear')
512            yb = bakInt(xdata)
513    if 'difC' in parmDict:
514        ff = 1.
515    else:       
516        try:
517            wave = parmDict[pfx+'Lam']
518        except KeyError:
519            wave = parmDict[pfx+'Lam1']
520        q = 4.0*np.pi*npsind(xdata/2.0)/wave
521        SQ = (q/(4.*np.pi))**2
522        FF = G2elem.GetFormFactorCoeff('Si')[0]
523        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
524    iD = 0       
525    while True:
526        try:
527            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
528            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
529            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
530            yb += ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
531            iD += 1       
532        except KeyError:
533            break
534    iD = 0
535    while True:
536        try:
537            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
538            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
539            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
540            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
541            shl = 0.002
542            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
543            iBeg = np.searchsorted(xdata,pkP-fmin)
544            iFin = np.searchsorted(xdata,pkP+fmax)
545            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
546            iD += 1       
547        except KeyError:
548            break
549        except ValueError:
550            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
551            break       
552    return yb
553   
554def getBackgroundDerv(pfx,parmDict,bakType,xdata):
555    'needs a doc string'
556    nBak = 0
557    while True:
558        key = pfx+'Back:'+str(nBak)
559        if key in parmDict:
560            nBak += 1
561        else:
562            break
563    dydb = np.zeros(shape=(nBak,len(xdata)))
564    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
565    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
566    cw = np.diff(xdata)
567    cw = np.append(cw,cw[-1])
568
569    if bakType in ['chebyschev','cosine']:
570        for iBak in range(nBak):   
571            if bakType == 'chebyschev':
572                dydb[iBak] = (xdata-xdata[0])**iBak
573            elif bakType == 'cosine':
574                dydb[iBak] = npcosd(xdata*iBak)
575    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
576        if nBak == 1:
577            dydb[0] = np.ones_like(xdata)
578        elif nBak == 2:
579            dX = xdata[-1]-xdata[0]
580            T2 = (xdata-xdata[0])/dX
581            T1 = 1.0-T2
582            dydb = [T1,T2]
583        else:
584            if bakType == 'lin interpolate':
585                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
586            elif bakType == 'inv interpolate':
587                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
588            elif bakType == 'log interpolate':
589                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
590            bakPos[0] = xdata[0]
591            bakPos[-1] = xdata[-1]
592            for i,pos in enumerate(bakPos):
593                if i == 0:
594                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
595                elif i == len(bakPos)-1:
596                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
597                else:
598                    dydb[i] = np.where(xdata>bakPos[i],
599                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
600                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
601    if 'difC' in parmDict:
602        ff = 1.
603    else:
604        try:
605            wave = parmDict[pfx+'Lam']
606        except KeyError:
607            wave = parmDict[pfx+'Lam1']
608        q = 4.0*np.pi*npsind(xdata/2.0)/wave
609        SQ = (q/(4*np.pi))**2
610        FF = G2elem.GetFormFactorCoeff('Si')[0]
611        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
612    iD = 0       
613    while True:
614        try:
615            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
616            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
617            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
618            sqr = np.sin(q*dbR)/(q*dbR)
619            cqr = np.cos(q*dbR)
620            temp = np.exp(-dbU*q**2)
621            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
622            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
623            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
624            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
625        except KeyError:
626            break
627    iD = 0
628    while True:
629        try:
630            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
631            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
632            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
633            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
634            shl = 0.002
635            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
636            iBeg = np.searchsorted(xdata,pkP-fmin)
637            iFin = np.searchsorted(xdata,pkP+fmax)
638            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
639            dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
640            dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
641            dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
642            dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
643            iD += 1       
644        except KeyError:
645            break
646        except ValueError:
647            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
648            break       
649    return dydb,dyddb,dydpk
650
651#use old fortran routine
652def getFCJVoigt3(pos,sig,gam,shl,xdata):
653    'needs a doc string'
654   
655    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
656#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
657    Df /= np.sum(Df)
658    return Df
659
660def getdFCJVoigt3(pos,sig,gam,shl,xdata):
661    'needs a doc string'
662   
663    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
664#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
665    sumDf = np.sum(Df)
666    return Df,dFdp,dFds,dFdg,dFdsh
667
668def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
669    'needs a doc string'
670    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
671    Df /= np.sum(Df)
672    return Df 
673   
674def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
675    'needs a doc string'
676    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
677    sumDf = np.sum(Df)
678    return Df,dFdp,dFda,dFdb,dFds,dFdg   
679
680def ellipseSize(H,Sij,GB):
681    'needs a doc string'
682    HX = np.inner(H.T,GB)
683    lenHX = np.sqrt(np.sum(HX**2))
684    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
685    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
686    lenR = np.sqrt(np.sum(R**2))
687    return lenR
688
689def ellipseSizeDerv(H,Sij,GB):
690    'needs a doc string'
691    lenR = ellipseSize(H,Sij,GB)
692    delt = 0.001
693    dRdS = np.zeros(6)
694    for i in range(6):
695        dSij = Sij[:]
696        dSij[i] += delt
697        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
698    return lenR,dRdS
699
700def getHKLpeak(dmin,SGData,A):
701    'needs a doc string'
702    HKL = G2lat.GenHLaue(dmin,SGData,A)       
703    HKLs = []
704    for h,k,l,d in HKL:
705        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
706        if not ext:
707            HKLs.append([h,k,l,d,-1])
708    return HKLs
709
710def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
711    'needs a doc string'
712   
713    yb = getBackground('',parmDict,bakType,xdata)
714    yc = np.zeros_like(yb)
715    cw = np.diff(xdata)
716    cw = np.append(cw,cw[-1])
717    if 'C' in dataType:
718        U = parmDict['U']
719        V = parmDict['V']
720        W = parmDict['W']
721        X = parmDict['X']
722        Y = parmDict['Y']
723        shl = max(parmDict['SH/L'],0.002)
724        Ka2 = False
725        if 'Lam1' in parmDict.keys():
726            Ka2 = True
727            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
728            kRatio = parmDict['I(L2)/I(L1)']
729        iPeak = 0
730        while True:
731            try:
732                pos = parmDict['pos'+str(iPeak)]
733                theta = (pos-parmDict['Zero'])/2.0
734                intens = parmDict['int'+str(iPeak)]
735                sigName = 'sig'+str(iPeak)
736                if sigName in varyList:
737                    sig = parmDict[sigName]
738                else:
739                    sig = U*tand(theta)**2+V*tand(theta)+W
740                sig = max(sig,0.001)          #avoid neg sigma
741                gamName = 'gam'+str(iPeak)
742                if gamName in varyList:
743                    gam = parmDict[gamName]
744                else:
745                    gam = X/cosd(theta)+Y*tand(theta)
746                gam = max(gam,0.001)             #avoid neg gamma
747                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
748                iBeg = np.searchsorted(xdata,pos-fmin)
749                iFin = np.searchsorted(xdata,pos+fmin)
750                if not iBeg+iFin:       #peak below low limit
751                    iPeak += 1
752                    continue
753                elif not iBeg-iFin:     #peak above high limit
754                    return yb+yc
755                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
756                if Ka2:
757                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
758                    iBeg = np.searchsorted(xdata,pos2-fmin)
759                    iFin = np.searchsorted(xdata,pos2+fmin)
760                    if iBeg-iFin:
761                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
762                iPeak += 1
763            except KeyError:        #no more peaks to process
764                return yb+yc
765    else:
766        Pdabc = parmDict['Pdabc']
767        difC = parmDict['difC']
768        alp0 = parmDict['alpha']
769        bet0 = parmDict['beta-0']
770        bet1 = parmDict['beta-1']
771        sig0 = parmDict['sig-0']
772        sig1 = parmDict['sig-1']
773        X = parmDict['X']
774        Y = parmDict['Y']
775        iPeak = 0
776        while True:
777            try:
778                pos = parmDict['pos'+str(iPeak)]               
779                tof = pos-parmDict['Zero']
780                dsp = tof/difC
781                intens = parmDict['int'+str(iPeak)]
782                alpName = 'alp'+str(iPeak)
783                if alpName in varyList:
784                    alp = parmDict[alpName]
785                else:
786                    if len(Pdabc):
787                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
788                    else:
789                        alp = alp0/dsp
790                betName = 'bet'+str(iPeak)
791                if betName in varyList:
792                    bet = parmDict[betName]
793                else:
794                    if len(Pdabc):
795                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
796                    else:
797                        bet = bet0+bet1/dsp**4
798                sigName = 'sig'+str(iPeak)
799                if sigName in varyList:
800                    sig = parmDict[sigName]
801                else:
802                    sig = sig0+sig1*dsp**2
803                gamName = 'gam'+str(iPeak)
804                if gamName in varyList:
805                    gam = parmDict[gamName]
806                else:
807                    gam = X*dsp**2+Y*dsp
808                gam = max(gam,0.001)             #avoid neg gamma
809                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
810                iBeg = np.searchsorted(xdata,pos-fmin)
811                iFin = np.searchsorted(xdata,pos+fmax)
812                lenX = len(xdata)
813                if not iBeg:
814                    iFin = np.searchsorted(xdata,pos+fmax)
815                elif iBeg == lenX:
816                    iFin = iBeg
817                else:
818                    iFin = np.searchsorted(xdata,pos+fmax)
819                if not iBeg+iFin:       #peak below low limit
820                    iPeak += 1
821                    continue
822                elif not iBeg-iFin:     #peak above high limit
823                    return yb+yc
824                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
825                iPeak += 1
826            except KeyError:        #no more peaks to process
827                return yb+yc
828           
829def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
830    'needs a doc string'
831# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
832    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
833    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
834    if 'Back:0' in varyList:            #background derivs are in front if present
835        dMdv[0:len(dMdb)] = dMdb
836    names = ['DebyeA','DebyeR','DebyeU']
837    for name in varyList:
838        if 'Debye' in name:
839            parm,id = name.split(':')
840            ip = names.index(parm)
841            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
842    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
843    for name in varyList:
844        if 'BkPk' in name:
845            parm,id = name.split(':')
846            ip = names.index(parm)
847            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
848    cw = np.diff(xdata)
849    cw = np.append(cw,cw[-1])
850    if 'C' in dataType:
851        U = parmDict['U']
852        V = parmDict['V']
853        W = parmDict['W']
854        X = parmDict['X']
855        Y = parmDict['Y']
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 = U*tanth**2+V*tanth+W
875                    dsdU = tanth**2
876                    dsdV = tanth
877                    dsdW = 1.0
878                sig = max(sig,0.001)          #avoid neg sigma
879                gamName = 'gam'+str(iPeak)
880                if gamName in varyList:
881                    gam = parmDict[gamName]
882                else:
883                    gam = X/costh+Y*tanth
884                    dgdX = 1.0/costh
885                    dgdY = tanth
886                gam = max(gam,0.001)             #avoid neg gamma
887                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
888                iBeg = np.searchsorted(xdata,pos-fmin)
889                iFin = np.searchsorted(xdata,pos+fmin)
890                if not iBeg+iFin:       #peak below low limit
891                    iPeak += 1
892                    continue
893                elif not iBeg-iFin:     #peak above high limit
894                    break
895                dMdpk = np.zeros(shape=(6,len(xdata)))
896                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
897                for i in range(1,5):
898                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
899                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
900                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
901                if Ka2:
902                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
903                    iBeg = np.searchsorted(xdata,pos2-fmin)
904                    iFin = np.searchsorted(xdata,pos2+fmin)
905                    if iBeg-iFin:
906                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
907                        for i in range(1,5):
908                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
909                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
910                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
911                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
912                for parmName in ['pos','int','sig','gam']:
913                    try:
914                        idx = varyList.index(parmName+str(iPeak))
915                        dMdv[idx] = dervDict[parmName]
916                    except ValueError:
917                        pass
918                if 'U' in varyList:
919                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
920                if 'V' in varyList:
921                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
922                if 'W' in varyList:
923                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
924                if 'X' in varyList:
925                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
926                if 'Y' in varyList:
927                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
928                if 'SH/L' in varyList:
929                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
930                if 'I(L2)/I(L1)' in varyList:
931                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
932                iPeak += 1
933            except KeyError:        #no more peaks to process
934                break
935    else:
936        Pdabc = parmDict['Pdabc']
937        difC = parmDict['difC']
938        alp0 = parmDict['alpha']
939        bet0 = parmDict['beta-0']
940        bet1 = parmDict['beta-1']
941        sig0 = parmDict['sig-0']
942        sig1 = parmDict['sig-1']
943        X = parmDict['X']
944        Y = parmDict['Y']
945        iPeak = 0
946        while True:
947            try:
948                pos = parmDict['pos'+str(iPeak)]               
949                tof = pos-parmDict['Zero']
950                dsp = tof/difC
951                intens = parmDict['int'+str(iPeak)]
952                alpName = 'alp'+str(iPeak)
953                if alpName in varyList:
954                    alp = parmDict[alpName]
955                else:
956                    if len(Pdabc):
957                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
958                    else:
959                        alp = alp0/dsp
960                    dada0 = 1./dsp
961                betName = 'bet'+str(iPeak)
962                if betName in varyList:
963                    bet = parmDict[betName]
964                else:
965                    if len(Pdabc):
966                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
967                    else:
968                        bet = bet0+bet1/dsp**4
969                    dbdb0 = 1.
970                    dbdb1 = 1/dsp**4
971                sigName = 'sig'+str(iPeak)
972                if sigName in varyList:
973                    sig = parmDict[sigName]
974                else:
975                    sig = sig0+sig1*dsp**2
976                    dsds0 = 1.
977                    dsds1 = dsp**2
978                gamName = 'gam'+str(iPeak)
979                if gamName in varyList:
980                    gam = parmDict[gamName]
981                else:
982                    gam = X*dsp**2+Y*dsp
983                    dsdX = dsp**2
984                    dsdY = dsp
985                gam = max(gam,0.001)             #avoid neg gamma
986                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
987                iBeg = np.searchsorted(xdata,pos-fmin)
988                lenX = len(xdata)
989                if not iBeg:
990                    iFin = np.searchsorted(xdata,pos+fmax)
991                elif iBeg == lenX:
992                    iFin = iBeg
993                else:
994                    iFin = np.searchsorted(xdata,pos+fmax)
995                if not iBeg+iFin:       #peak below low limit
996                    iPeak += 1
997                    continue
998                elif not iBeg-iFin:     #peak above high limit
999                    break
1000                dMdpk = np.zeros(shape=(7,len(xdata)))
1001                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1002                for i in range(1,6):
1003                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1004                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1005                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1006                for parmName in ['pos','int','alp','bet','sig','gam']:
1007                    try:
1008                        idx = varyList.index(parmName+str(iPeak))
1009                        dMdv[idx] = dervDict[parmName]
1010                    except ValueError:
1011                        pass
1012                if 'alpha' in varyList:
1013                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1014                if 'beta-0' in varyList:
1015                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1016                if 'beta-1' in varyList:
1017                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1018                if 'sig-0' in varyList:
1019                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1020                if 'sig-1' in varyList:
1021                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1022                if 'X' in varyList:
1023                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1024                if 'Y' in varyList:
1025                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1026                iPeak += 1
1027            except KeyError:        #no more peaks to process
1028                break
1029    return dMdv
1030       
1031def Dict2Values(parmdict, varylist):
1032    '''Use before call to leastsq to setup list of values for the parameters
1033    in parmdict, as selected by key in varylist'''
1034    return [parmdict[key] for key in varylist] 
1035   
1036def Values2Dict(parmdict, varylist, values):
1037    ''' Use after call to leastsq to update the parameter dictionary with
1038    values corresponding to keys in varylist'''
1039    parmdict.update(zip(varylist,values))
1040   
1041def SetBackgroundParms(Background):
1042    'needs a doc string'
1043    if len(Background) == 1:            # fix up old backgrounds
1044        BackGround.append({'nDebye':0,'debyeTerms':[]})
1045    bakType,bakFlag = Background[0][:2]
1046    backVals = Background[0][3:]
1047    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1048    Debye = Background[1]           #also has background peaks stuff
1049    backDict = dict(zip(backNames,backVals))
1050    backVary = []
1051    if bakFlag:
1052        backVary = backNames
1053
1054    backDict['nDebye'] = Debye['nDebye']
1055    debyeDict = {}
1056    debyeList = []
1057    for i in range(Debye['nDebye']):
1058        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1059        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1060        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1061    debyeVary = []
1062    for item in debyeList:
1063        if item[1]:
1064            debyeVary.append(item[0])
1065    backDict.update(debyeDict)
1066    backVary += debyeVary
1067
1068    backDict['nPeaks'] = Debye['nPeaks']
1069    peaksDict = {}
1070    peaksList = []
1071    for i in range(Debye['nPeaks']):
1072        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1073        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1074        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1075    peaksVary = []
1076    for item in peaksList:
1077        if item[1]:
1078            peaksVary.append(item[0])
1079    backDict.update(peaksDict)
1080    backVary += peaksVary   
1081    return bakType,backDict,backVary
1082           
1083def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1084    'needs a doc string'
1085   
1086       
1087    def GetBackgroundParms(parmList,Background):
1088        iBak = 0
1089        while True:
1090            try:
1091                bakName = 'Back:'+str(iBak)
1092                Background[0][iBak+3] = parmList[bakName]
1093                iBak += 1
1094            except KeyError:
1095                break
1096        iDb = 0
1097        while True:
1098            names = ['DebyeA:','DebyeR:','DebyeU:']
1099            try:
1100                for i,name in enumerate(names):
1101                    val = parmList[name+str(iDb)]
1102                    Background[1]['debyeTerms'][iDb][2*i] = val
1103                iDb += 1
1104            except KeyError:
1105                break
1106        iDb = 0
1107        while True:
1108            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1109            try:
1110                for i,name in enumerate(names):
1111                    val = parmList[name+str(iDb)]
1112                    Background[1]['peaksList'][iDb][2*i] = val
1113                iDb += 1
1114            except KeyError:
1115                break
1116               
1117    def BackgroundPrint(Background,sigDict):
1118        if Background[0][1]:
1119            print 'Background coefficients for',Background[0][0],'function'
1120            ptfmt = "%12.5f"
1121            ptstr =  'values:'
1122            sigstr = 'esds  :'
1123            for i,back in enumerate(Background[0][3:]):
1124                ptstr += ptfmt % (back)
1125                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1126            print ptstr
1127            print sigstr
1128        else:
1129            print 'Background not refined'
1130        if Background[1]['nDebye']:
1131            parms = ['DebyeA','DebyeR','DebyeU']
1132            print 'Debye diffuse scattering coefficients'
1133            ptfmt = "%12.5f"
1134            names =   'names :'
1135            ptstr =  'values:'
1136            sigstr = 'esds  :'
1137            for item in sigDict:
1138                if 'Debye' in item:
1139                    names += '%12s'%(item)
1140                    sigstr += ptfmt%(sigDict[item])
1141                    parm,id = item.split(':')
1142                    ip = parms.index(parm)
1143                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1144            print names
1145            print ptstr
1146            print sigstr
1147        if Background[1]['nPeaks']:
1148            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1149            print 'Peaks in background coefficients'
1150            ptfmt = "%12.5f"
1151            names =   'names :'
1152            ptstr =  'values:'
1153            sigstr = 'esds  :'
1154            for item in sigDict:
1155                if 'BkPk' in item:
1156                    names += '%12s'%(item)
1157                    sigstr += ptfmt%(sigDict[item])
1158                    parm,id = item.split(':')
1159                    ip = parms.index(parm)
1160                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1161            print names
1162            print ptstr
1163            print sigstr
1164                           
1165    def SetInstParms(Inst):
1166        dataType = Inst['Type'][0]
1167        insVary = []
1168        insNames = []
1169        insVals = []
1170        for parm in Inst:
1171            insNames.append(parm)
1172            insVals.append(Inst[parm][1])
1173            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1174                'beta-0','beta-1','sig-0','sig-1',] and Inst[parm][2]:
1175                    insVary.append(parm)
1176        instDict = dict(zip(insNames,insVals))
1177        instDict['X'] = max(instDict['X'],0.01)
1178        instDict['Y'] = max(instDict['Y'],0.01)
1179        if 'SH/L' in instDict:
1180            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1181        return dataType,instDict,insVary
1182       
1183    def GetInstParms(parmDict,Inst,varyList,Peaks):
1184        for name in Inst:
1185            Inst[name][1] = parmDict[name]
1186        iPeak = 0
1187        while True:
1188            try:
1189                sigName = 'sig'+str(iPeak)
1190                pos = parmDict['pos'+str(iPeak)]
1191                if sigName not in varyList:
1192                    if 'C' in Inst['Type'][0]:
1193                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1194                    else:
1195                        dsp = pos/Inst['difC'][1]
1196                        parmDict[sigName] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1197                gamName = 'gam'+str(iPeak)
1198                if gamName not in varyList:
1199                    if 'C' in Inst['Type'][0]:
1200                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1201                    else:
1202                        dsp = pos/Inst['difC'][1]
1203                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1204                iPeak += 1
1205            except KeyError:
1206                break
1207       
1208    def InstPrint(Inst,sigDict):
1209        print 'Instrument Parameters:'
1210        ptfmt = "%12.6f"
1211        ptlbls = 'names :'
1212        ptstr =  'values:'
1213        sigstr = 'esds  :'
1214        for parm in Inst:
1215            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1216                'beta-0','beta-1','sig-0','sig-1',]:
1217                ptlbls += "%s" % (parm.center(12))
1218                ptstr += ptfmt % (Inst[parm][1])
1219                if parm in sigDict:
1220                    sigstr += ptfmt % (sigDict[parm])
1221                else:
1222                    sigstr += 12*' '
1223        print ptlbls
1224        print ptstr
1225        print sigstr
1226
1227    def SetPeaksParms(dataType,Peaks):
1228        peakNames = []
1229        peakVary = []
1230        peakVals = []
1231        if 'C' in dataType:
1232            names = ['pos','int','sig','gam']
1233        else:
1234            names = ['pos','int','alp','bet','sig','gam']
1235        for i,peak in enumerate(Peaks):
1236            for j,name in enumerate(names):
1237                peakVals.append(peak[2*j])
1238                parName = name+str(i)
1239                peakNames.append(parName)
1240                if peak[2*j+1]:
1241                    peakVary.append(parName)
1242        return dict(zip(peakNames,peakVals)),peakVary
1243               
1244    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1245        if 'C' in Inst['Type'][0]:
1246            names = ['pos','int','sig','gam']
1247        else:
1248            names = ['pos','int','alp','bet','sig','gam']
1249        for i,peak in enumerate(Peaks):
1250            pos = parmDict['pos'+str(i)]
1251            if 'difC' in Inst:
1252                dsp = pos/Inst['difC'][1]
1253            for j in range(len(names)):
1254                parName = names[j]+str(i)
1255                if parName in varyList:
1256                    peak[2*j] = parmDict[parName]
1257                elif 'alpha' in parName:
1258                    peak[2*j] = parmDict['alpha']/dsp
1259                elif 'beta' in parName:
1260                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
1261                elif 'sig' in parName:
1262                    if 'C' in Inst['Type'][0]:
1263                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1264                    else:
1265                        peak[2*j] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1266                elif 'gam' in parName:
1267                    if 'C' in Inst['Type'][0]:
1268                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1269                    else:
1270                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1271                       
1272    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1273        print 'Peak coefficients:'
1274        if 'C' in dataType:
1275            names = ['pos','int','sig','gam']
1276        else:
1277            names = ['pos','int','alp','bet','sig','gam']           
1278        head = 13*' '
1279        for name in names:
1280            head += name.center(10)+'esd'.center(10)
1281        print head
1282        if 'C' in dataType:
1283            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1284        else:
1285            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1286        for i,peak in enumerate(Peaks):
1287            ptstr =  ':'
1288            for j in range(len(names)):
1289                name = names[j]
1290                parName = name+str(i)
1291                ptstr += ptfmt[name] % (parmDict[parName])
1292                if parName in varyList:
1293#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1294                    ptstr += ptfmt[name] % (sigDict[parName])
1295                else:
1296#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1297                    ptstr += 10*' '
1298            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1299               
1300    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1301        parmdict.update(zip(varylist,values))
1302        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1303           
1304    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1305        parmdict.update(zip(varylist,values))
1306        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1307        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1308        if dlg:
1309            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1310            if not GoOn:
1311                return -M           #abort!!
1312        return M
1313       
1314    if controls:
1315        Ftol = controls['min dM/M']
1316        derivType = controls['deriv type']
1317    else:
1318        Ftol = 0.0001
1319        derivType = 'analytic'
1320    if oneCycle:
1321        Ftol = 1.0
1322    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1323    yc *= 0.                            #set calcd ones to zero
1324    yb *= 0.
1325    yd *= 0.
1326    xBeg = np.searchsorted(x,Limits[0])
1327    xFin = np.searchsorted(x,Limits[1])
1328    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1329    dataType,insDict,insVary = SetInstParms(Inst)
1330    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1331    parmDict = {}
1332    parmDict.update(bakDict)
1333    parmDict.update(insDict)
1334    parmDict.update(peakDict)
1335    parmDict['Pdabc'] = []      #dummy Pdabc
1336    parmDict.update(Inst2)      #put in real one if there
1337    varyList = bakVary+insVary+peakVary
1338    while True:
1339        begin = time.time()
1340        values =  np.array(Dict2Values(parmDict, varyList))
1341        if FitPgm == 'LSQ':
1342            try:
1343                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1344                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1345                ncyc = int(result[2]['nfev']/2)
1346            finally:
1347                dlg.Destroy()
1348            runtime = time.time()-begin   
1349            chisq = np.sum(result[2]['fvec']**2)
1350            Values2Dict(parmDict, varyList, result[0])
1351            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1352            GOF = chisq/(xFin-xBeg-len(varyList))
1353            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1354            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1355            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1356            try:
1357                sig = np.sqrt(np.diag(result[1])*GOF)
1358                if np.any(np.isnan(sig)):
1359                    print '*** Least squares aborted - some invalid esds possible ***'
1360                break                   #refinement succeeded - finish up!
1361            except ValueError:          #result[1] is None on singular matrix
1362                print '**** Refinement failed - singular matrix ****'
1363                Ipvt = result[2]['ipvt']
1364                for i,ipvt in enumerate(Ipvt):
1365                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1366                        print 'Removing parameter: ',varyList[ipvt-1]
1367                        del(varyList[ipvt-1])
1368                        break
1369        elif FitPgm == 'BFGS':
1370            print 'Other program here'
1371            return
1372       
1373    sigDict = dict(zip(varyList,sig))
1374    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1375    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1376    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1377    GetBackgroundParms(parmDict,Background)
1378    BackgroundPrint(Background,sigDict)
1379    GetInstParms(parmDict,Inst,varyList,Peaks)
1380    InstPrint(Inst,sigDict)
1381    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1382    PeaksPrint(dataType,parmDict,sigDict,varyList)
1383
1384def calcIncident(Iparm,xdata):
1385    'needs a doc string'
1386
1387    def IfunAdv(Iparm,xdata):
1388        Itype = Iparm['Itype']
1389        Icoef = Iparm['Icoeff']
1390        DYI = np.ones((12,xdata.shape[0]))
1391        YI = np.ones_like(xdata)*Icoef[0]
1392       
1393        x = xdata/1000.                 #expressions are in ms
1394        if Itype == 'Exponential':
1395            for i in range(1,10,2):
1396                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1397                YI += Icoef[i]*Eterm
1398                DYI[i] *= Eterm
1399                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1400        elif 'Maxwell'in Itype:
1401            Eterm = np.exp(-Icoef[2]/x**2)
1402            DYI[1] = Eterm/x**5
1403            DYI[2] = -Icoef[1]*DYI[1]/x**2
1404            YI += (Icoef[1]*Eterm/x**5)
1405            if 'Exponential' in Itype:
1406                for i in range(3,12,2):
1407                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1408                    YI += Icoef[i]*Eterm
1409                    DYI[i] *= Eterm
1410                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1411            else:   #Chebyschev
1412                T = (2./x)-1.
1413                Ccof = np.ones((12,xdata.shape[0]))
1414                Ccof[1] = T
1415                for i in range(2,12):
1416                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1417                for i in range(1,10):
1418                    YI += Ccof[i]*Icoef[i+2]
1419                    DYI[i+2] =Ccof[i]
1420        return YI,DYI
1421       
1422    Iesd = np.array(Iparm['Iesd'])
1423    Icovar = Iparm['Icovar']
1424    YI,DYI = IfunAdv(Iparm,xdata)
1425    YI = np.where(YI>0,YI,1.)
1426    WYI = np.zeros_like(xdata)
1427    vcov = np.zeros((12,12))
1428    k = 0
1429    for i in range(12):
1430        for j in range(i,12):
1431            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1432            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1433            k += 1
1434    M = np.inner(vcov,DYI.T)
1435    WYI = np.sum(M*DYI,axis=0)
1436    WYI = np.where(WYI>0.,WYI,0.)
1437    return YI,WYI
1438   
1439#testing data
1440NeedTestData = True
1441def TestData():
1442    'needs a doc string'
1443#    global NeedTestData
1444    NeedTestData = False
1445    global bakType
1446    bakType = 'chebyschev'
1447    global xdata
1448    xdata = np.linspace(4.0,40.0,36000)
1449    global parmDict0
1450    parmDict0 = {
1451        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1452        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1453        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1454        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1455        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1456        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1457        }
1458    global parmDict1
1459    parmDict1 = {
1460        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1461        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1462        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1463        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1464        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1465        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1466        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1467        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1468        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1469        }
1470    global parmDict2
1471    parmDict2 = {
1472        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1473        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1474        'Back0':5.,'Back1':-0.02,'Back2':.004,
1475#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1476        }
1477    global varyList
1478    varyList = []
1479
1480def test0():
1481    if NeedTestData: TestData()
1482    msg = 'test '
1483    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1484    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1485    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1486    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1487    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1488    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1489   
1490def test1():
1491    if NeedTestData: TestData()
1492    time0 = time.time()
1493    for i in range(100):
1494        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1495    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1496   
1497def test2(name,delt):
1498    if NeedTestData: TestData()
1499    varyList = [name,]
1500    xdata = np.linspace(5.6,5.8,400)
1501    hplot = plotter.add('derivatives test for '+name).gca()
1502    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1503    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1504    parmDict2[name] += delt
1505    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1506    hplot.plot(xdata,(y1-y0)/delt,'r+')
1507   
1508def test3(name,delt):
1509    if NeedTestData: TestData()
1510    names = ['pos','sig','gam','shl']
1511    idx = names.index(name)
1512    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1513    xdata = np.linspace(5.6,5.8,800)
1514    dx = xdata[1]-xdata[0]
1515    hplot = plotter.add('derivatives test for '+name).gca()
1516    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1517    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1518    myDict[name] += delt
1519    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1520    hplot.plot(xdata,(y1-y0)/delt,'r+')
1521
1522if __name__ == '__main__':
1523    import GSASIItestplot as plot
1524    global plotter
1525    plotter = plot.PlotNotebook()
1526#    test0()
1527#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1528    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1529        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1530        test2(name,shft)
1531    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1532        test3(name,shft)
1533    print "OK"
1534    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.