source: trunk/GSASIIpwd.py @ 1073

Last change on this file since 1073 was 1073, checked in by vondreele, 8 years ago

change chebyschev background fxn to match what is in old GSAS. Much better fits now with it.

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