source: trunk/GSASIIpwd.py @ 1893

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

fix background fitting by interpolation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 69.6 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-17 16:47:40 +0000 (Wed, 17 Jun 2015) $
10# $Author: toby $
11# $Revision: 1893 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1893 2015-06-17 16:47:40Z 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: 1893 $")
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        print 'Background coefficients for',Background[0][0],'function'
1380        ptfmt = "%12.5f"
1381        ptstr =  'value: '
1382        sigstr = 'esd  : '
1383        for i,back in enumerate(Background[0][3:]):
1384            ptstr += ptfmt % (back)
1385            if Background[0][1]:
1386                sigstr += ptfmt % (sigDict['Back;'+str(i)])
1387            if len(ptstr) > 75:
1388                print ptstr
1389                if Background[0][1]: print sigstr
1390                ptstr =  'value: '
1391                sigstr = 'esd  : '
1392        if len(ptstr) > 8:
1393            print ptstr
1394            if Background[0][1]: print sigstr
1395
1396        if Background[1]['nDebye']:
1397            parms = ['DebyeA;','DebyeR;','DebyeU;']
1398            print 'Debye diffuse scattering coefficients'
1399            ptfmt = "%12.5f"
1400            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1401            for term in range(Background[1]['nDebye']):
1402                line = ' term %d'%(term)
1403                for ip,name in enumerate(parms):
1404                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1405                    if name+str(term) in sigDict:
1406                        line += ptfmt%(sigDict[name+str(term)])
1407                print line
1408        if Background[1]['nPeaks']:
1409            print 'Coefficients for Background Peaks'
1410            ptfmt = "%15.3f"
1411            for j,pl in enumerate(Background[1]['peaksList']):
1412                names =  'peak %3d:'%(j+1)
1413                ptstr =  'values  :'
1414                sigstr = 'esds    :'
1415                for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']):
1416                    val = pl[2*i]
1417                    prm = lbl+";"+str(j)
1418                    names += '%15s'%(prm)
1419                    ptstr += ptfmt%(val)
1420                    if prm in sigDict:
1421                        sigstr += ptfmt%(sigDict[prm])
1422                    else:
1423                        sigstr += " "*15
1424                print names
1425                print ptstr
1426                print sigstr
1427                           
1428    def SetInstParms(Inst):
1429        dataType = Inst['Type'][0]
1430        insVary = []
1431        insNames = []
1432        insVals = []
1433        for parm in Inst:
1434            insNames.append(parm)
1435            insVals.append(Inst[parm][1])
1436            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1437                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1438                    insVary.append(parm)
1439        instDict = dict(zip(insNames,insVals))
1440        instDict['X'] = max(instDict['X'],0.01)
1441        instDict['Y'] = max(instDict['Y'],0.01)
1442        if 'SH/L' in instDict:
1443            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1444        return dataType,instDict,insVary
1445       
1446    def GetInstParms(parmDict,Inst,varyList,Peaks):
1447        for name in Inst:
1448            Inst[name][1] = parmDict[name]
1449        iPeak = 0
1450        while True:
1451            try:
1452                sigName = 'sig'+str(iPeak)
1453                pos = parmDict['pos'+str(iPeak)]
1454                if sigName not in varyList:
1455                    if 'C' in Inst['Type'][0]:
1456                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1457                    else:
1458                        dsp = G2lat.Pos2dsp(Inst,pos)
1459                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1460                gamName = 'gam'+str(iPeak)
1461                if gamName not in varyList:
1462                    if 'C' in Inst['Type'][0]:
1463                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1464                    else:
1465                        dsp = G2lat.Pos2dsp(Inst,pos)
1466                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1467                iPeak += 1
1468            except KeyError:
1469                break
1470       
1471    def InstPrint(Inst,sigDict):
1472        print 'Instrument Parameters:'
1473        ptfmt = "%12.6f"
1474        ptlbls = 'names :'
1475        ptstr =  'values:'
1476        sigstr = 'esds  :'
1477        for parm in Inst:
1478            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1479                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1480                ptlbls += "%s" % (parm.center(12))
1481                ptstr += ptfmt % (Inst[parm][1])
1482                if parm in sigDict:
1483                    sigstr += ptfmt % (sigDict[parm])
1484                else:
1485                    sigstr += 12*' '
1486        print ptlbls
1487        print ptstr
1488        print sigstr
1489
1490    def SetPeaksParms(dataType,Peaks):
1491        peakNames = []
1492        peakVary = []
1493        peakVals = []
1494        if 'C' in dataType:
1495            names = ['pos','int','sig','gam']
1496        else:
1497            names = ['pos','int','alp','bet','sig','gam']
1498        for i,peak in enumerate(Peaks):
1499            for j,name in enumerate(names):
1500                peakVals.append(peak[2*j])
1501                parName = name+str(i)
1502                peakNames.append(parName)
1503                if peak[2*j+1]:
1504                    peakVary.append(parName)
1505        return dict(zip(peakNames,peakVals)),peakVary
1506               
1507    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1508        if 'C' in Inst['Type'][0]:
1509            names = ['pos','int','sig','gam']
1510        else:   #'T'
1511            names = ['pos','int','alp','bet','sig','gam']
1512        for i,peak in enumerate(Peaks):
1513            pos = parmDict['pos'+str(i)]
1514            if 'difC' in Inst:
1515                dsp = pos/Inst['difC'][1]
1516            for j in range(len(names)):
1517                parName = names[j]+str(i)
1518                if parName in varyList:
1519                    peak[2*j] = parmDict[parName]
1520                elif 'alpha' in parName:
1521                    peak[2*j] = parmDict['alpha']/dsp
1522                elif 'beta' in parName:
1523                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1524                elif 'sig' in parName:
1525                    if 'C' in Inst['Type'][0]:
1526                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1527                    else:
1528                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1529                elif 'gam' in parName:
1530                    if 'C' in Inst['Type'][0]:
1531                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1532                    else:
1533                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1534                       
1535    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1536        print 'Peak coefficients:'
1537        if 'C' in dataType:
1538            names = ['pos','int','sig','gam']
1539        else:   #'T'
1540            names = ['pos','int','alp','bet','sig','gam']           
1541        head = 13*' '
1542        for name in names:
1543            if name in ['alp','bet']:
1544                head += name.center(8)+'esd'.center(8)
1545            else:
1546                head += name.center(10)+'esd'.center(10)
1547        print head
1548        if 'C' in dataType:
1549            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1550        else:
1551            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1552        for i,peak in enumerate(Peaks):
1553            ptstr =  ':'
1554            for j in range(len(names)):
1555                name = names[j]
1556                parName = name+str(i)
1557                ptstr += ptfmt[name] % (parmDict[parName])
1558                if parName in varyList:
1559                    ptstr += ptfmt[name] % (sigDict[parName])
1560                else:
1561                    if name in ['alp','bet']:
1562                        ptstr += 8*' '
1563                    else:
1564                        ptstr += 10*' '
1565            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1566               
1567    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1568        parmdict.update(zip(varylist,values))
1569        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1570           
1571    def errPeakProfile(values,xdata,ydata,weights,dataType,parmdict,varylist,bakType,dlg):       
1572        parmdict.update(zip(varylist,values))
1573        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1574        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1575        if dlg:
1576            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1577            if not GoOn:
1578                return -M           #abort!!
1579        return M
1580       
1581    if controls:
1582        Ftol = controls['min dM/M']
1583        derivType = controls['deriv type']
1584    else:
1585        Ftol = 0.0001
1586        derivType = 'analytic'
1587    if oneCycle:
1588        Ftol = 1.0
1589    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1590    yc *= 0.                            #set calcd ones to zero
1591    yb *= 0.
1592    yd *= 0.
1593    xBeg = np.searchsorted(x,Limits[0])
1594    xFin = np.searchsorted(x,Limits[1])
1595    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1596    dataType,insDict,insVary = SetInstParms(Inst)
1597    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1598    parmDict = {}
1599    parmDict.update(bakDict)
1600    parmDict.update(insDict)
1601    parmDict.update(peakDict)
1602    parmDict['Pdabc'] = []      #dummy Pdabc
1603    parmDict.update(Inst2)      #put in real one if there
1604    if prevVaryList:
1605        varyList = prevVaryList[:]
1606    else:
1607        varyList = bakVary+insVary+peakVary
1608    fullvaryList = varyList[:]
1609    while True:
1610        begin = time.time()
1611        values =  np.array(Dict2Values(parmDict, varyList))
1612        Rvals = {}
1613        badVary = []
1614        if FitPgm == 'LSQ':
1615            try:
1616                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1617                   args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1618                ncyc = int(result[2]['nfev']/2)
1619            finally:
1620                if dlg: dlg.Destroy()
1621            runtime = time.time()-begin   
1622            chisq = np.sum(result[2]['fvec']**2)
1623            Values2Dict(parmDict, varyList, result[0])
1624            Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1625            Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1626            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1627            if ncyc:
1628                print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1629            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1630            sig = [0]*len(varyList)
1631            if len(varyList) == 0: break  # if nothing was refined
1632            try:
1633                sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1634                if np.any(np.isnan(sig)):
1635                    print '*** Least squares aborted - some invalid esds possible ***'
1636                break                   #refinement succeeded - finish up!
1637            except ValueError:          #result[1] is None on singular matrix
1638                print '**** Refinement failed - singular matrix ****'
1639                Ipvt = result[2]['ipvt']
1640                for i,ipvt in enumerate(Ipvt):
1641                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1642                        print 'Removing parameter: ',varyList[ipvt-1]
1643                        badVary.append(varyList[ipvt-1])
1644                        del(varyList[ipvt-1])
1645                        break
1646                else: # nothing removed
1647                    break
1648        elif FitPgm == 'BFGS':
1649            print 'Other program here'
1650            return {}
1651       
1652    sigDict = dict(zip(varyList,sig))
1653    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1654    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1655    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1656    GetBackgroundParms(parmDict,Background)
1657    if bakVary: BackgroundPrint(Background,sigDict)
1658    GetInstParms(parmDict,Inst,varyList,Peaks)
1659    if insVary: InstPrint(Inst,sigDict)
1660    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1661    if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList)
1662    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1663
1664def calcIncident(Iparm,xdata):
1665    'needs a doc string'
1666
1667    def IfunAdv(Iparm,xdata):
1668        Itype = Iparm['Itype']
1669        Icoef = Iparm['Icoeff']
1670        DYI = np.ones((12,xdata.shape[0]))
1671        YI = np.ones_like(xdata)*Icoef[0]
1672       
1673        x = xdata/1000.                 #expressions are in ms
1674        if Itype == 'Exponential':
1675            for i in [1,3,5,7,9]:
1676                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1677                YI += Icoef[i]*Eterm
1678                DYI[i] *= Eterm
1679                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1680        elif 'Maxwell'in Itype:
1681            Eterm = np.exp(-Icoef[2]/x**2)
1682            DYI[1] = Eterm/x**5
1683            DYI[2] = -Icoef[1]*DYI[1]/x**2
1684            YI += (Icoef[1]*Eterm/x**5)
1685            if 'Exponential' in Itype:
1686                for i in range(3,12,2):
1687                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1688                    YI += Icoef[i]*Eterm
1689                    DYI[i] *= Eterm
1690                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1691            else:   #Chebyschev
1692                T = (2./x)-1.
1693                Ccof = np.ones((12,xdata.shape[0]))
1694                Ccof[1] = T
1695                for i in range(2,12):
1696                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1697                for i in range(1,10):
1698                    YI += Ccof[i]*Icoef[i+2]
1699                    DYI[i+2] =Ccof[i]
1700        return YI,DYI
1701       
1702    Iesd = np.array(Iparm['Iesd'])
1703    Icovar = Iparm['Icovar']
1704    YI,DYI = IfunAdv(Iparm,xdata)
1705    YI = np.where(YI>0,YI,1.)
1706    WYI = np.zeros_like(xdata)
1707    vcov = np.zeros((12,12))
1708    k = 0
1709    for i in range(12):
1710        for j in range(i,12):
1711            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1712            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1713            k += 1
1714    M = np.inner(vcov,DYI.T)
1715    WYI = np.sum(M*DYI,axis=0)
1716    WYI = np.where(WYI>0.,WYI,0.)
1717    return YI,WYI
1718   
1719#testing data
1720NeedTestData = True
1721def TestData():
1722    'needs a doc string'
1723#    global NeedTestData
1724    NeedTestData = False
1725    global bakType
1726    bakType = 'chebyschev'
1727    global xdata
1728    xdata = np.linspace(4.0,40.0,36000)
1729    global parmDict0
1730    parmDict0 = {
1731        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1732        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1733        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1734        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1735        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1736        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1737        }
1738    global parmDict1
1739    parmDict1 = {
1740        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1741        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1742        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1743        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1744        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1745        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1746        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1747        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1748        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1749        }
1750    global parmDict2
1751    parmDict2 = {
1752        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1753        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1754        'Back0':5.,'Back1':-0.02,'Back2':.004,
1755#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1756        }
1757    global varyList
1758    varyList = []
1759
1760def test0():
1761    if NeedTestData: TestData()
1762    msg = 'test '
1763    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1764    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
1765    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1766    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1767    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
1768    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1769   
1770def test1():
1771    if NeedTestData: TestData()
1772    time0 = time.time()
1773    for i in range(100):
1774        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1775    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1776   
1777def test2(name,delt):
1778    if NeedTestData: TestData()
1779    varyList = [name,]
1780    xdata = np.linspace(5.6,5.8,400)
1781    hplot = plotter.add('derivatives test for '+name).gca()
1782    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1783    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1784    parmDict2[name] += delt
1785    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1786    hplot.plot(xdata,(y1-y0)/delt,'r+')
1787   
1788def test3(name,delt):
1789    if NeedTestData: TestData()
1790    names = ['pos','sig','gam','shl']
1791    idx = names.index(name)
1792    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1793    xdata = np.linspace(5.6,5.8,800)
1794    dx = xdata[1]-xdata[0]
1795    hplot = plotter.add('derivatives test for '+name).gca()
1796    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1797    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1798    myDict[name] += delt
1799    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1800    hplot.plot(xdata,(y1-y0)/delt,'r+')
1801
1802if __name__ == '__main__':
1803    import GSASIItestplot as plot
1804    global plotter
1805    plotter = plot.PlotNotebook()
1806#    test0()
1807#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1808    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1809        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1810        test2(name,shft)
1811    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1812        test3(name,shft)
1813    print "OK"
1814    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.