source: trunk/GSASIIpwd.py @ 1889

Last change on this file since 1889 was 1889, checked in by toby, 7 years ago

Add fixed background points; animate dragging points and lines; add fixed bkg fit routine (bkg peaks buggy\!)

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