source: trunk/GSASIIpwd.py @ 1874

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

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