source: trunk/GSASIIpwd.py @ 1682

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

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