source: trunk/GSASIIpwd.py @ 798

Last change on this file since 798 was 798, checked in by vondreele, 13 years ago

more peak fitting work including TOF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.1 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2012-11-14 22:37:26 +0000 (Wed, 14 Nov 2012) $
6# $Author: vondreele $
7# $Revision: 798 $
8# $URL: trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 798 2012-11-14 22:37:26Z vondreele $
10########### SVN repository information ###################
11import sys
12import math
13import time
14
15import numpy as np
16import scipy as sp
17import numpy.linalg as nl
18from numpy.fft import ifft, fft, fftshift
19import scipy.interpolate as si
20import scipy.stats as st
21import scipy.optimize as so
22
23import GSASIIpath
24GSASIIpath.SetVersionNumber("$Revision: 798 $")
25import GSASIIlattice as G2lat
26import GSASIIspc as G2spc
27import GSASIIElem as G2elem
28import GSASIIgrid as G2gd
29import GSASIIIO as G2IO
30import GSASIImath as G2mth
31import pypowder as pyd
32
33# trig functions in degrees
34sind = lambda x: math.sin(x*math.pi/180.)
35asind = lambda x: 180.*math.asin(x)/math.pi
36tand = lambda x: math.tan(x*math.pi/180.)
37atand = lambda x: 180.*math.atan(x)/math.pi
38atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
39cosd = lambda x: math.cos(x*math.pi/180.)
40acosd = lambda x: 180.*math.acos(x)/math.pi
41rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
42#numpy versions
43npsind = lambda x: np.sin(x*np.pi/180.)
44npasind = lambda x: 180.*np.arcsin(x)/math.pi
45npcosd = lambda x: np.cos(x*math.pi/180.)
46npacosd = lambda x: 180.*np.arccos(x)/math.pi
47nptand = lambda x: np.tan(x*math.pi/180.)
48npatand = lambda x: 180.*np.arctan(x)/np.pi
49npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
51npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
52   
53#GSASII pdf calculation routines
54       
55def Transmission(Geometry,Abs,Diam):
56#Calculate sample transmission
57#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate'
58#   Abs: absorption coeff in cm-1
59#   Diam: sample thickness/diameter in mm
60    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
61        MuR = Abs*Diam/20.0
62        if MuR <= 3.0:
63            T0 = 16/(3.*math.pi)
64            T1 = -0.045780
65            T2 = -0.02489
66            T3 = 0.003045
67            T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
68            if T < -20.:
69                return 2.06e-9
70            else:
71                return math.exp(T)
72        else:
73            T1 = 1.433902
74            T2 = 0.013869+0.337894
75            T3 = 1.933433+1.163198
76            T4 = 0.044365-0.04259
77            T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
78            return T/100.
79    elif 'plate' in Geometry:
80        MuR = Abs*Diam/10.
81        return math.exp(-MuR)
82    elif 'Bragg' in Geometry:
83        return 0.0
84
85def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0):
86#Calculate sample absorption
87#   Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate'
88#   MuR: absorption coeff * sample thickness/2 or radius
89#   Tth: 2-theta scattering angle - can be numpy array
90#   Phi: flat plate tilt angle - future
91#   Psi: flat plate tilt axis - future
92    Sth2 = npsind(Tth/2.0)**2
93    Cth2 = 1.-Sth2
94    if 'Cylinder' in Geometry:      #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample
95        if MuR < 3.0:
96            T0 = 16.0/(3*np.pi)
97            T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \
98                0.109561*np.sqrt(Sth2)-26.04556
99            T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \
100                1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3
101            T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2
102            Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4
103            return np.exp(Trns)
104        else:
105            T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \
106                10.02088*Sth2**3-3.36778*Sth2**4
107            T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \
108                (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039
109            T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \
110                0.13576*np.sqrt(Sth2)+1.163198
111            T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202
112            Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4
113            return Trns/100.
114    elif 'Bragg' in Geometry:
115        return 1.0
116    elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam
117        # and only defined for 2theta < 90
118        MuT = 2.*MuR
119        T1 = np.exp(-MuT)
120        T2 = np.exp(-MuT/npcosd(Tth))
121        Tb = MuT-MuT/npcosd(Tth)
122        return (T2-T1)/Tb
123    elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector
124        MuT = 2.*MuR
125        cth = npcosd(Tth/2.0)
126        return np.exp(-MuT/cth)/cth
127       
128def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0):
129    dA = 0.001
130    AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi)
131    if MuR:
132        AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi)
133        return (AbsP-AbsM)/(2.0*dA)
134    else:
135        return (AbsP-1.)/dA
136       
137def Polarization(Pola,Tth,Azm=0.0):
138    """   Calculate angle dependent x-ray polarization correction (not scaled correctly!)
139    input:
140        Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized
141        Azm: azimuthal angle e.g. 0.0 in plane of polarization
142        Tth: 2-theta scattering angle - can be numpy array
143            which (if either) of these is "right"?
144    return:
145        pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \
146            (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
147        dpdPola: derivative needed for least squares
148    """
149    pola = ((1.0-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+   \
150        (1.0-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2
151    dpdPola = -npsind(Tth)**2*(npsind(Azm)**2-npcosd(Azm)**2)
152    return pola,dpdPola
153   
154def Oblique(ObCoeff,Tth):
155    if ObCoeff:
156        return (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth)))
157    else:
158        return 1.0
159               
160def Ruland(RulCoff,wave,Q,Compton):
161    C = 2.9978e8
162    D = 1.5e-3
163    hmc = 0.024262734687
164    sinth2 = (Q*wave/(4.0*np.pi))**2
165    dlam = (wave**2)*Compton*Q/C
166    dlam_c = 2.0*hmc*sinth2-D*wave**2
167    return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2))
168   
169def LorchWeight(Q):
170    return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1]))
171           
172def GetAsfMean(ElList,Sthl2):
173#   Calculate various scattering factor terms for PDF calcs
174#   ElList: element dictionary contains scattering factor coefficients, etc.
175#   Sthl2: numpy array of sin theta/lambda squared values
176#   returns: mean(f^2), mean(f)^2, mean(compton)
177    sumNoAtoms = 0.0
178    FF = np.zeros_like(Sthl2)
179    FF2 = np.zeros_like(Sthl2)
180    CF = np.zeros_like(Sthl2)
181    for El in ElList:
182        sumNoAtoms += ElList[El]['FormulaNo']
183    for El in ElList:
184        el = ElList[El]
185        ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2
186        cf = G2elem.ComptonFac(el,Sthl2)
187        FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms
188        FF2 += ff2*el['FormulaNo']/sumNoAtoms
189        CF += cf*el['FormulaNo']/sumNoAtoms
190    return FF2,FF**2,CF
191   
192def GetNumDensity(ElList,Vol):
193    sumNoAtoms = 0.0
194    for El in ElList:
195        sumNoAtoms += ElList[El]['FormulaNo']
196    return sumNoAtoms/Vol
197           
198def CalcPDF(data,inst,xydata):
199    auxPlot = []
200    import copy
201    import scipy.fftpack as ft
202    #subtract backgrounds - if any
203    xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
204    if data['Sample Bkg.']['Name']:
205        xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1]+
206            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
207    if data['Container']['Name']:
208        xycontainer = (xydata['Container'][1][1]+data['Container']['Add'])*data['Container']['Mult']
209        if data['Container Bkg.']['Name']:
210            xycontainer += (xydata['Container Bkg.'][1][1]+
211                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
212        xydata['IofQ'][1][1] += xycontainer
213    #get element data & absorption coeff.
214    ElList = data['ElList']
215    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
216    #Apply angle dependent corrections
217    Tth = xydata['Sample'][1][0]
218    dt = (Tth[1]-Tth[0])
219    MuR = Abs*data['Diam']/20.0
220    xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
221    xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
222    if data['DetType'] == 'Image plate':
223        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
224    XY = xydata['IofQ'][1]   
225    #convert to Q
226    hc = 12.397639
227    wave = G2mth.getWave(inst)
228    keV = hc/wave
229    minQ = npT2q(Tth[0],wave)
230    maxQ = npT2q(Tth[-1],wave)   
231    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
232    dq = Qpoints[1]-Qpoints[0]
233    XY[0] = npT2q(XY[0],wave)   
234#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
235    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
236    Qdata = T(Qpoints)
237   
238    qLimits = data['QScaleLim']
239    minQ = np.searchsorted(Qpoints,qLimits[0])
240    maxQ = np.searchsorted(Qpoints,qLimits[1])
241    newdata = []
242    xydata['IofQ'][1][0] = Qpoints
243    xydata['IofQ'][1][1] = Qdata
244    for item in xydata['IofQ'][1]:
245        newdata.append(item[:maxQ])
246    xydata['IofQ'][1] = newdata
247   
248
249    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
250    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
251    Q = xydata['SofQ'][1][0]
252    ruland = Ruland(data['Ruland'],wave,Q,CF)
253#    auxPlot.append([Q,ruland,'Ruland'])     
254    CF *= ruland
255#    auxPlot.append([Q,CF,'CF-Corr'])
256    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
257    xydata['SofQ'][1][1] *= scale
258    xydata['SofQ'][1][1] -= CF
259    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
260    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
261    xydata['SofQ'][1][1] *= scale
262   
263    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
264    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
265    if data['Lorch']:
266        xydata['FofQ'][1][1] *= LorchWeight(Q)
267   
268    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
269    nR = len(xydata['GofR'][1][1])
270    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
271    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
272   
273       
274    return auxPlot
275       
276#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
277
278def factorize(num):
279    ''' Provide prime number factors for integer num
280    Returns dictionary of prime factors (keys) & power for each (data)
281    '''
282    factors = {}
283    orig = num
284
285    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
286    i, sqi = 2, 4
287    while sqi <= num:
288        while not num%i:
289            num /= i
290            factors[i] = factors.get(i, 0) + 1
291
292        sqi += 2*i + 1
293        i += 1
294
295    if num != 1 and num != orig:
296        factors[num] = factors.get(num, 0) + 1
297
298    if factors:
299        return factors
300    else:
301        return {num:1}          #a prime number!
302           
303def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
304    ''' Provide list of optimal data sizes for FFT calculations
305    Input:
306        nmin: minimum data size >= 1
307        nmax: maximum data size > nmin
308        thresh: maximum prime factor allowed
309    Returns:
310        list of data sizes where the maximum prime factor is < thresh
311    ''' 
312    plist = []
313    nmin = max(1,nmin)
314    nmax = max(nmin+1,nmax)
315    for p in range(nmin,nmax):
316        if max(factorize(p).keys()) < thresh:
317            plist.append(p)
318    return plist
319
320np.seterr(divide='ignore')
321
322# Normal distribution
323
324# loc = mu, scale = std
325_norm_pdf_C = 1./math.sqrt(2*math.pi)
326class norm_gen(st.rv_continuous):
327       
328    def pdf(self,x,*args,**kwds):
329        loc,scale=kwds['loc'],kwds['scale']
330        x = (x-loc)/scale
331        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
332       
333norm = norm_gen(name='norm',longname='A normal',extradoc="""
334
335Normal distribution
336
337The location (loc) keyword specifies the mean.
338The scale (scale) keyword specifies the standard deviation.
339
340normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
341""")
342
343## Cauchy
344
345# median = loc
346
347class cauchy_gen(st.rv_continuous):
348
349    def pdf(self,x,*args,**kwds):
350        loc,scale=kwds['loc'],kwds['scale']
351        x = (x-loc)/scale
352        return 1.0/np.pi/(1.0+x*x) / scale
353       
354cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
355
356Cauchy distribution
357
358cauchy.pdf(x) = 1/(pi*(1+x**2))
359
360This is the t distribution with one degree of freedom.
361""")
362   
363   
364#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
365
366
367class fcjde_gen(st.rv_continuous):
368    """
369    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
370    Ref: J. Appl. Cryst. (1994) 27, 892-900.
371    Parameters
372    -----------------------------------------
373    x: array -1 to 1
374    t: 2-theta position of peak
375    s: sum(S/L,H/L); S: sample height, H: detector opening,
376        L: sample to detector opening distance
377    dx: 2-theta step size in deg
378    Result for fcj.pdf
379    -----------------------------------------
380    T = x*dx+t
381    s = S/L+H/L
382    if x < 0:
383        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|   
384    if x >= 0:
385        fcj.pdf = 0   
386    """
387    def _pdf(self,x,t,s,dx):
388        T = dx*x+t
389        ax2 = abs(npcosd(T))
390        ax = ax2**2
391        bx = npcosd(t)**2
392        bx = np.where(ax>bx,bx,ax)
393        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
394        fx = np.where(fx > 0.,fx,0.0)
395        return fx
396             
397    def pdf(self,x,*args,**kwds):
398        loc=kwds['loc']
399        return self._pdf(x-loc,*args)
400       
401fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
402               
403def getWidthsCW(pos,sig,gam,shl):
404    widths = [np.sqrt(sig)/100.,gam/200.]
405    fwhm = 2.355*widths[0]+2.*widths[1]
406    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
407    fmax = 15.0*fwhm
408    if pos > 90:
409        fmin,fmax = [fmax,fmin]         
410    return widths,fmin,fmax
411   
412def getWidthsTOF(pos,alp,bet,sig,gam):
413    lnf = 3.3      # =log(0.001/2)
414    widths = [np.sqrt(sig),gam]
415    fwhm = 2.355*widths[0]+2.*widths[1]
416    fmin = 10.*fwhm*(1.+1./alp)   
417    fmax = 10.*fwhm*(1.+1./bet)
418    return widths,fmin,fmax
419   
420def getFWHM(TTh,Inst):
421    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
422    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
423    gamFW = lambda s,g: math.exp(math.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.)
424    s = sig(TTh/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
425    g = gam(TTh/2.,Inst['X'][1],Inst['Y'][1])*100.
426    return gamFW(g,s)   
427               
428def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
429    DX = xdata[1]-xdata[0]
430    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
431    x = np.linspace(pos-fmin,pos+fmin,256)
432    dx = x[1]-x[0]
433    Norm = norm.pdf(x,loc=pos,scale=widths[0])
434    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
435    arg = [pos,shl/57.2958,dx,]
436    FCJ = fcjde.pdf(x,*arg,loc=pos)
437    if len(np.nonzero(FCJ)[0])>5:
438        z = np.column_stack([Norm,Cauchy,FCJ]).T
439        Z = fft(z)
440        Df = ifft(Z.prod(axis=0)).real
441    else:
442        z = np.column_stack([Norm,Cauchy]).T
443        Z = fft(z)
444        Df = fftshift(ifft(Z.prod(axis=0))).real
445    Df /= np.sum(Df)
446    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
447    return intens*Df(xdata)*DX/dx
448
449def getBackground(pfx,parmDict,bakType,xdata):
450    yb = np.zeros_like(xdata)
451    nBak = 0
452    while True:
453        key = pfx+'Back:'+str(nBak)
454        if key in parmDict:
455            nBak += 1
456        else:
457            break
458    if bakType in ['chebyschev','cosine']:
459        for iBak in range(nBak):   
460            key = pfx+'Back:'+str(iBak)
461            if bakType == 'chebyschev':
462                yb += parmDict[key]*(xdata-xdata[0])**iBak
463            elif bakType == 'cosine':
464                yb += parmDict[key]*npcosd(xdata*iBak)
465    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
466        if nBak == 1:
467            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
468        elif nBak == 2:
469            dX = xdata[-1]-xdata[0]
470            T2 = (xdata-xdata[0])/dX
471            T1 = 1.0-T2
472            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
473        else:
474            if bakType == 'lin interpolate':
475                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
476            elif bakType == 'inv interpolate':
477                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
478            elif bakType == 'log interpolate':
479                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
480            bakPos[0] = xdata[0]
481            bakPos[-1] = xdata[-1]
482            bakVals = np.zeros(nBak)
483            for i in range(nBak):
484                bakVals[i] = parmDict[pfx+'Back:'+str(i)]
485            bakInt = si.interp1d(bakPos,bakVals,'linear')
486            yb = bakInt(xdata)
487    if 'difC' in parmDict:
488        ff = 1.
489    else:       
490        try:
491            wave = parmDict[pfx+'Lam']
492        except KeyError:
493            wave = parmDict[pfx+'Lam1']
494        q = 4.0*np.pi*npsind(xdata/2.0)/wave
495        SQ = (q/(4.*np.pi))**2
496        FF = G2elem.GetFormFactorCoeff('Si')[0]
497        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
498    iD = 0       
499    while True:
500        try:
501            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
502            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
503            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
504            yb += ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
505            iD += 1       
506        except KeyError:
507            break
508    iD = 0
509    while True:
510        try:
511            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
512            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
513            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
514            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
515            shl = 0.002
516            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
517            iBeg = np.searchsorted(xdata,pkP-fmin)
518            iFin = np.searchsorted(xdata,pkP+fmax)
519            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
520            iD += 1       
521        except KeyError:
522            break       
523    return yb
524   
525def getBackgroundDerv(pfx,parmDict,bakType,xdata):
526    nBak = 0
527    while True:
528        key = pfx+'Back:'+str(nBak)
529        if key in parmDict:
530            nBak += 1
531        else:
532            break
533    dydb = np.zeros(shape=(nBak,len(xdata)))
534    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
535    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
536    dx = xdata[1]-xdata[0]
537
538    if bakType in ['chebyschev','cosine']:
539        for iBak in range(nBak):   
540            if bakType == 'chebyschev':
541                dydb[iBak] = (xdata-xdata[0])**iBak
542            elif bakType == 'cosine':
543                dydb[iBak] = npcosd(xdata*iBak)
544    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
545        if nBak == 1:
546            dydb[0] = np.ones_like(xdata)
547        elif nBak == 2:
548            dX = xdata[-1]-xdata[0]
549            T2 = (xdata-xdata[0])/dX
550            T1 = 1.0-T2
551            dydb = [T1,T2]
552        else:
553            if bakType == 'lin interpolate':
554                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
555            elif bakType == 'inv interpolate':
556                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
557            elif bakType == 'log interpolate':
558                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
559            bakPos[0] = xdata[0]
560            bakPos[-1] = xdata[-1]
561            dx = bakPos[1]-bakPos[0]
562            for i,pos in enumerate(bakPos):
563                if i == 0:
564                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
565                elif i == len(bakPos)-1:
566                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
567                else:
568                    dydb[i] = np.where(xdata>bakPos[i],
569                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
570                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
571    if 'difC' in parmDict:
572        ff = 1.
573    else:
574        try:
575            wave = parmDict[pfx+'Lam']
576        except KeyError:
577            wave = parmDict[pfx+'Lam1']
578        q = 4.0*np.pi*npsind(xdata/2.0)/wave
579        SQ = (q/(4*np.pi))**2
580        FF = G2elem.GetFormFactorCoeff('Si')[0]
581        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
582    iD = 0       
583    while True:
584        try:
585            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
586            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
587            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
588            sqr = np.sin(q*dbR)/(q*dbR)
589            cqr = np.cos(q*dbR)
590            temp = np.exp(-dbU*q**2)
591            dyddb[3*iD] = ff*sqr*temp/(np.pi*dx)
592            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*dx)
593            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*dx)
594            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
595        except KeyError:
596            break
597    iD = 0
598    while True:
599        try:
600            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
601            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
602            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
603            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
604            shl = 0.002
605            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
606            iBeg = np.searchsorted(xdata,pkP-fmin)
607            iFin = np.searchsorted(xdata,pkP+fmax)
608            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
609            dydpk[4*iD][iBeg:iFin] += 100.*dx*pkI*dFdp
610            dydpk[4*iD+1][iBeg:iFin] += 100.*dx*Df
611            dydpk[4*iD+2][iBeg:iFin] += 100.*dx*pkI*dFds
612            dydpk[4*iD+3][iBeg:iFin] += 100.*dx*pkI*dFdg
613            iD += 1       
614        except KeyError:
615            break
616    return dydb,dyddb,dydpk
617
618#use old fortran routine
619def getFCJVoigt3(pos,sig,gam,shl,xdata):
620   
621    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
622#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
623    Df /= np.sum(Df)
624    return Df
625
626def getdFCJVoigt3(pos,sig,gam,shl,xdata):
627   
628    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
629#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
630    sumDf = np.sum(Df)
631    return Df,dFdp,dFds,dFdg,dFdsh
632
633def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
634    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
635    Df /= np.sum(Df)
636    return Df 
637   
638def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
639    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
640    sumDf = np.sum(Df)
641    return Df,dFdp,dFda,dFdb,dFds,dFdg   
642
643def ellipseSize(H,Sij,GB):
644    HX = np.inner(H.T,GB)
645    lenHX = np.sqrt(np.sum(HX**2))
646    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
647    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
648    lenR = np.sqrt(np.sum(R**2))
649    return lenR
650
651def ellipseSizeDerv(H,Sij,GB):
652    lenR = ellipseSize(H,Sij,GB)
653    delt = 0.001
654    dRdS = np.zeros(6)
655    for i in range(6):
656        dSij = Sij[:]
657        dSij[i] += delt
658        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
659    return lenR,dRdS
660
661def getHKLpeak(dmin,SGData,A):
662    HKL = G2lat.GenHLaue(dmin,SGData,A)       
663    HKLs = []
664    for h,k,l,d in HKL:
665        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
666        if not ext:
667            HKLs.append([h,k,l,d,-1])
668    return HKLs
669
670def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
671   
672    yb = getBackground('',parmDict,bakType,xdata)
673    yc = np.zeros_like(yb)
674    if 'C' in dataType:
675        dx = xdata[1]-xdata[0]
676        U = parmDict['U']
677        V = parmDict['V']
678        W = parmDict['W']
679        X = parmDict['X']
680        Y = parmDict['Y']
681        shl = max(parmDict['SH/L'],0.002)
682        Ka2 = False
683        if 'Lam1' in parmDict.keys():
684            Ka2 = True
685            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
686            kRatio = parmDict['I(L2)/I(L1)']
687        iPeak = 0
688        while True:
689            try:
690                pos = parmDict['pos'+str(iPeak)]
691                theta = (pos-parmDict['Zero'])/2.0
692                intens = parmDict['int'+str(iPeak)]
693                sigName = 'sig'+str(iPeak)
694                if sigName in varyList:
695                    sig = parmDict[sigName]
696                else:
697                    sig = U*tand(theta)**2+V*tand(theta)+W
698                sig = max(sig,0.001)          #avoid neg sigma
699                gamName = 'gam'+str(iPeak)
700                if gamName in varyList:
701                    gam = parmDict[gamName]
702                else:
703                    gam = X/cosd(theta)+Y*tand(theta)
704                gam = max(gam,0.001)             #avoid neg gamma
705                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
706                iBeg = np.searchsorted(xdata,pos-fmin)
707                lenX = len(xdata)
708                if not iBeg:
709                    iFin = np.searchsorted(xdata,pos+fmin)
710                elif iBeg == lenX:
711                    iFin = iBeg
712                else:
713                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
714                if not iBeg+iFin:       #peak below low limit
715                    iPeak += 1
716                    continue
717                elif not iBeg-iFin:     #peak above high limit
718                    return yb+yc
719                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
720                if Ka2:
721                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
722                    kdelt = int((pos2-pos)/dx)               
723                    iBeg = min(lenX,iBeg+kdelt)
724                    iFin = min(lenX,iFin+kdelt)
725                    if iBeg-iFin:
726                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
727                iPeak += 1
728            except KeyError:        #no more peaks to process
729                return yb+yc
730    else:
731        Pdabc = parmDict['Pdabc']
732        difC = parmDict['difC']
733        alp0 = parmDict['alpha']
734        bet0 = parmDict['beta-0']
735        bet1 = parmDict['beta-1']
736        sig0 = parmDict['sig-0']
737        sig1 = parmDict['sig-1']
738        X = parmDict['X']
739        Y = parmDict['Y']
740        iPeak = 0
741        while True:
742            try:
743                pos = parmDict['pos'+str(iPeak)]               
744                tof = pos-parmDict['Zero']
745                dsp = tof/difC
746                intens = parmDict['int'+str(iPeak)]
747                alpName = 'alp'+str(iPeak)
748                if alpName in varyList:
749                    alp = parmDict[alpName]
750                else:
751                    if len(Pdabc):
752                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
753                    else:
754                        alp = alp0/dsp
755                betName = 'bet'+str(iPeak)
756                if betName in varyList:
757                    bet = parmDict[betName]
758                else:
759                    if len(Pdabc):
760                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
761                    else:
762                        bet = bet0+bet1/dsp**4
763                sigName = 'sig'+str(iPeak)
764                if sigName in varyList:
765                    sig = parmDict[sigName]
766                else:
767                    sig = sig0+sig1*dsp**2
768                gamName = 'gam'+str(iPeak)
769                if gamName in varyList:
770                    gam = parmDict[gamName]
771                else:
772                    gam = X*dsp**2+Y*dsp
773                gam = max(gam,0.001)             #avoid neg gamma
774                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
775                iBeg = np.searchsorted(xdata,pos-fmin)
776                iFin = np.searchsorted(xdata,pos+fmax)
777                lenX = len(xdata)
778                if not iBeg:
779                    iFin = np.searchsorted(xdata,pos+fmax)
780                elif iBeg == lenX:
781                    iFin = iBeg
782                else:
783                    iFin = np.searchsorted(xdata,pos+fmax)
784                if not iBeg+iFin:       #peak below low limit
785                    iPeak += 1
786                    continue
787                elif not iBeg-iFin:     #peak above high limit
788                    return yb+yc
789                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
790                iPeak += 1
791            except KeyError:        #no more peaks to process
792                return yb+yc
793           
794def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
795# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
796    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
797    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
798    if 'Back:0' in varyList:            #background derivs are in front if present
799        dMdv[0:len(dMdb)] = dMdb
800    names = ['DebyeA','DebyeR','DebyeU']
801    for name in varyList:
802        if 'Debye' in name:
803            parm,id = name.split(':')
804            ip = names.index(parm)
805            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
806    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
807    for name in varyList:
808        if 'BkPk' in name:
809            parm,id = name.split(':')
810            ip = names.index(parm)
811            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
812    if 'C' in dataType:
813        dx = xdata[1]-xdata[0]
814        U = parmDict['U']
815        V = parmDict['V']
816        W = parmDict['W']
817        X = parmDict['X']
818        Y = parmDict['Y']
819        shl = max(parmDict['SH/L'],0.002)
820        Ka2 = False
821        if 'Lam1' in parmDict.keys():
822            Ka2 = True
823            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
824            kRatio = parmDict['I(L2)/I(L1)']
825        iPeak = 0
826        while True:
827            try:
828                pos = parmDict['pos'+str(iPeak)]
829                theta = (pos-parmDict['Zero'])/2.0
830                intens = parmDict['int'+str(iPeak)]
831                sigName = 'sig'+str(iPeak)
832                tanth = tand(theta)
833                costh = cosd(theta)
834                if sigName in varyList:
835                    sig = parmDict[sigName]
836                else:
837                    sig = U*tanth**2+V*tanth+W
838                    dsdU = tanth**2
839                    dsdV = tanth
840                    dsdW = 1.0
841                sig = max(sig,0.001)          #avoid neg sigma
842                gamName = 'gam'+str(iPeak)
843                if gamName in varyList:
844                    gam = parmDict[gamName]
845                else:
846                    gam = X/costh+Y*tanth
847                    dgdX = 1.0/costh
848                    dgdY = tanth
849                gam = max(gam,0.001)             #avoid neg gamma
850                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
851                iBeg = np.searchsorted(xdata,pos-fmin)
852                lenX = len(xdata)
853                if not iBeg:
854                    iFin = np.searchsorted(xdata,pos+fmin)
855                elif iBeg == lenX:
856                    iFin = iBeg
857                else:
858                    iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
859                if not iBeg+iFin:       #peak below low limit
860                    iPeak += 1
861                    continue
862                elif not iBeg-iFin:     #peak above high limit
863                    break
864                dMdpk = np.zeros(shape=(6,len(xdata)))
865                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
866                for i in range(1,5):
867                    dMdpk[i][iBeg:iFin] += 100.*dx*intens*dMdipk[i]
868                dMdpk[0][iBeg:iFin] += 100.*dx*dMdipk[0]
869                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
870                if Ka2:
871                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
872                    kdelt = int((pos2-pos)/dx)               
873                    iBeg = min(lenX,iBeg+kdelt)
874                    iFin = min(lenX,iFin+kdelt)
875                    if iBeg-iFin:
876                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
877                        for i in range(1,5):
878                            dMdpk[i][iBeg:iFin] += 100.*dx*intens*kRatio*dMdipk2[i]
879                        dMdpk[0][iBeg:iFin] += 100.*dx*kRatio*dMdipk2[0]
880                        dMdpk[5][iBeg:iFin] += 100.*dx*dMdipk2[0]
881                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
882                for parmName in ['pos','int','sig','gam']:
883                    try:
884                        idx = varyList.index(parmName+str(iPeak))
885                        dMdv[idx] = dervDict[parmName]
886                    except ValueError:
887                        pass
888                if 'U' in varyList:
889                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
890                if 'V' in varyList:
891                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
892                if 'W' in varyList:
893                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
894                if 'X' in varyList:
895                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
896                if 'Y' in varyList:
897                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
898                if 'SH/L' in varyList:
899                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
900                if 'I(L2)/I(L1)' in varyList:
901                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
902                iPeak += 1
903            except KeyError:        #no more peaks to process
904                break
905    else:
906        Pdabc = parmDict['Pdabc']
907        difC = parmDict['difC']
908        alp0 = parmDict['alpha']
909        bet0 = parmDict['beta-0']
910        bet1 = parmDict['beta-1']
911        sig0 = parmDict['sig-0']
912        sig1 = parmDict['sig-1']
913        X = parmDict['X']
914        Y = parmDict['Y']
915        iPeak = 0
916        while True:
917            try:
918                pos = parmDict['pos'+str(iPeak)]               
919                tof = pos-parmDict['Zero']
920                dsp = tof/difC
921                intens = parmDict['int'+str(iPeak)]
922                alpName = 'alp'+str(iPeak)
923                if alpName in varyList:
924                    alp = parmDict[alpName]
925                else:
926                    if len(Pdabc):
927                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
928                    else:
929                        alp = alp0/dsp
930                    dada0 = 1./dsp
931                betName = 'bet'+str(iPeak)
932                if betName in varyList:
933                    bet = parmDict[betName]
934                else:
935                    if len(Pdabc):
936                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
937                    else:
938                        bet = bet0+bet1/dsp**4
939                    dbdb0 = 1.
940                    dbdb1 = 1/dsp**4
941                sigName = 'sig'+str(iPeak)
942                if sigName in varyList:
943                    sig = parmDict[sigName]
944                else:
945                    sig = sig0+sig1*dsp**2
946                    dsds0 = 1.
947                    dsds1 = dsp**2
948                gamName = 'gam'+str(iPeak)
949                if gamName in varyList:
950                    gam = parmDict[gamName]
951                else:
952                    gam = X*dsp**2+Y*dsp
953                    dsdX = dsp**2
954                    dsdY = dsp
955                gam = max(gam,0.001)             #avoid neg gamma
956                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
957                iBeg = np.searchsorted(xdata,pos-fmin)
958                lenX = len(xdata)
959                if not iBeg:
960                    iFin = np.searchsorted(xdata,pos+fmax)
961                elif iBeg == lenX:
962                    iFin = iBeg
963                else:
964                    iFin = np.searchsorted(xdata,pos+fmax)
965                if not iBeg+iFin:       #peak below low limit
966                    iPeak += 1
967                    continue
968                elif not iBeg-iFin:     #peak above high limit
969                    break
970                dMdpk = np.zeros(shape=(7,len(xdata)))
971                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
972                for i in range(1,6):
973                    dMdpk[i][iBeg:iFin] += intens*dMdipk[i]
974                dMdpk[0][iBeg:iFin] += dMdipk[0]
975                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
976                for parmName in ['pos','int','alp','bet','sig','gam']:
977                    try:
978                        idx = varyList.index(parmName+str(iPeak))
979                        dMdv[idx] = dervDict[parmName]
980                    except ValueError:
981                        pass
982                if 'alpha' in varyList:
983                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
984                if 'beta-0' in varyList:
985                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
986                if 'beta-1' in varyList:
987                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
988                if 'sig-0' in varyList:
989                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
990                if 'sig-1' in varyList:
991                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
992                if 'X' in varyList:
993                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
994                if 'Y' in varyList:
995                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
996                iPeak += 1
997            except KeyError:        #no more peaks to process
998                break
999    return dMdv
1000       
1001def Dict2Values(parmdict, varylist):
1002    '''Use before call to leastsq to setup list of values for the parameters
1003    in parmdict, as selected by key in varylist'''
1004    return [parmdict[key] for key in varylist] 
1005   
1006def Values2Dict(parmdict, varylist, values):
1007    ''' Use after call to leastsq to update the parameter dictionary with
1008    values corresponding to keys in varylist'''
1009    parmdict.update(zip(varylist,values))
1010   
1011def SetBackgroundParms(Background):
1012    if len(Background) == 1:            # fix up old backgrounds
1013        BackGround.append({'nDebye':0,'debyeTerms':[]})
1014    bakType,bakFlag = Background[0][:2]
1015    backVals = Background[0][3:]
1016    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1017    Debye = Background[1]           #also has background peaks stuff
1018    backDict = dict(zip(backNames,backVals))
1019    backVary = []
1020    if bakFlag:
1021        backVary = backNames
1022
1023    backDict['nDebye'] = Debye['nDebye']
1024    debyeDict = {}
1025    debyeList = []
1026    for i in range(Debye['nDebye']):
1027        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1028        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1029        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1030    debyeVary = []
1031    for item in debyeList:
1032        if item[1]:
1033            debyeVary.append(item[0])
1034    backDict.update(debyeDict)
1035    backVary += debyeVary
1036
1037    backDict['nPeaks'] = Debye['nPeaks']
1038    peaksDict = {}
1039    peaksList = []
1040    for i in range(Debye['nPeaks']):
1041        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1042        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1043        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1044    peaksVary = []
1045    for item in peaksList:
1046        if item[1]:
1047            peaksVary.append(item[0])
1048    backDict.update(peaksDict)
1049    backVary += peaksVary   
1050    return bakType,backDict,backVary
1051           
1052def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1053   
1054       
1055    def GetBackgroundParms(parmList,Background):
1056        iBak = 0
1057        while True:
1058            try:
1059                bakName = 'Back:'+str(iBak)
1060                Background[0][iBak+3] = parmList[bakName]
1061                iBak += 1
1062            except KeyError:
1063                break
1064        iDb = 0
1065        while True:
1066            names = ['DebyeA:','DebyeR:','DebyeU:']
1067            try:
1068                for i,name in enumerate(names):
1069                    val = parmList[name+str(iDb)]
1070                    Background[1]['debyeTerms'][iDb][2*i] = val
1071                iDb += 1
1072            except KeyError:
1073                break
1074        iDb = 0
1075        while True:
1076            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1077            try:
1078                for i,name in enumerate(names):
1079                    val = parmList[name+str(iDb)]
1080                    Background[1]['peaksList'][iDb][2*i] = val
1081                iDb += 1
1082            except KeyError:
1083                break
1084               
1085    def BackgroundPrint(Background,sigDict):
1086        if Background[0][1]:
1087            print 'Background coefficients for',Background[0][0],'function'
1088            ptfmt = "%12.5f"
1089            ptstr =  'values:'
1090            sigstr = 'esds  :'
1091            for i,back in enumerate(Background[0][3:]):
1092                ptstr += ptfmt % (back)
1093                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1094            print ptstr
1095            print sigstr
1096        else:
1097            print 'Background not refined'
1098        if Background[1]['nDebye']:
1099            parms = ['DebyeA','DebyeR','DebyeU']
1100            print 'Debye diffuse scattering coefficients'
1101            ptfmt = "%12.5f"
1102            names =   'names :'
1103            ptstr =  'values:'
1104            sigstr = 'esds  :'
1105            for item in sigDict:
1106                if 'Debye' in item:
1107                    names += '%12s'%(item)
1108                    sigstr += ptfmt%(sigDict[item])
1109                    parm,id = item.split(':')
1110                    ip = parms.index(parm)
1111                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1112            print names
1113            print ptstr
1114            print sigstr
1115        if Background[1]['nPeaks']:
1116            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1117            print 'Peaks in background coefficients'
1118            ptfmt = "%12.5f"
1119            names =   'names :'
1120            ptstr =  'values:'
1121            sigstr = 'esds  :'
1122            for item in sigDict:
1123                if 'BkPk' in item:
1124                    names += '%12s'%(item)
1125                    sigstr += ptfmt%(sigDict[item])
1126                    parm,id = item.split(':')
1127                    ip = parms.index(parm)
1128                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1129            print names
1130            print ptstr
1131            print sigstr
1132                           
1133    def SetInstParms(Inst):
1134        dataType = Inst['Type'][0]
1135        insVary = []
1136        insNames = []
1137        insVals = []
1138        for parm in Inst:
1139            insNames.append(parm)
1140            insVals.append(Inst[parm][1])
1141            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1142                'beta-0','beta-1','sig-0','sig-1',] and Inst[parm][2]:
1143                    insVary.append(parm)
1144        instDict = dict(zip(insNames,insVals))
1145        instDict['X'] = max(instDict['X'],0.01)
1146        instDict['Y'] = max(instDict['Y'],0.01)
1147        if 'SH/L' in instDict:
1148            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1149        return dataType,instDict,insVary
1150       
1151    def GetInstParms(parmDict,Inst,varyList,Peaks):
1152        for name in Inst:
1153            Inst[name][1] = parmDict[name]
1154        iPeak = 0
1155        while True:
1156            try:
1157                sigName = 'sig'+str(iPeak)
1158                pos = parmDict['pos'+str(iPeak)]
1159                if sigName not in varyList:
1160                    if 'C' in Inst['Type'][0]:
1161                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1162                    else:
1163                        dsp = pos/Inst['difC'][1]
1164                        parmDict[sigName] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1165                gamName = 'gam'+str(iPeak)
1166                if gamName not in varyList:
1167                    if 'C' in Inst['Type'][0]:
1168                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1169                    else:
1170                        dsp = pos/Inst['difC'][1]
1171                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1172                iPeak += 1
1173            except KeyError:
1174                break
1175       
1176    def InstPrint(Inst,sigDict):
1177        print 'Instrument Parameters:'
1178        ptfmt = "%12.6f"
1179        ptlbls = 'names :'
1180        ptstr =  'values:'
1181        sigstr = 'esds  :'
1182        for parm in Inst:
1183            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1184                'beta-0','beta-1','sig-0','sig-1',]:
1185                ptlbls += "%s" % (parm.center(12))
1186                ptstr += ptfmt % (Inst[parm][1])
1187                if parm in sigDict:
1188                    sigstr += ptfmt % (sigDict[parm])
1189                else:
1190                    sigstr += 12*' '
1191        print ptlbls
1192        print ptstr
1193        print sigstr
1194
1195    def SetPeaksParms(dataType,Peaks):
1196        peakNames = []
1197        peakVary = []
1198        peakVals = []
1199        if 'C' in dataType:
1200            names = ['pos','int','sig','gam']
1201        else:
1202            names = ['pos','int','alp','bet','sig','gam']
1203        for i,peak in enumerate(Peaks):
1204            for j,name in enumerate(names):
1205                peakVals.append(peak[2*j])
1206                parName = name+str(i)
1207                peakNames.append(parName)
1208                if peak[2*j+1]:
1209                    peakVary.append(parName)
1210        return dict(zip(peakNames,peakVals)),peakVary
1211               
1212    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1213        if 'C' in Inst['Type'][0]:
1214            names = ['pos','int','sig','gam']
1215        else:
1216            names = ['pos','int','alp','bet','sig','gam']
1217        for i,peak in enumerate(Peaks):
1218            pos = parmDict['pos'+str(i)]
1219            if 'difC' in Inst:
1220                dsp = pos/Inst['difC'][1]
1221            for j in range(len(names)):
1222                parName = names[j]+str(i)
1223                if parName in varyList:
1224                    peak[2*j] = parmDict[parName]
1225                elif 'alpha' in parName:
1226                    peak[2*j] = parmDict['alpha']/dsp
1227                elif 'beta' in parName:
1228                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
1229                elif 'sig' in parName:
1230                    if 'C' in Inst['Type'][0]:
1231                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1232                    else:
1233                        peak[2*j] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1234                elif 'gam' in parName:
1235                    if 'C' in Inst['Type'][0]:
1236                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1237                    else:
1238                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1239                       
1240    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1241        print 'Peak coefficients:'
1242        if 'C' in dataType:
1243            names = ['pos','int','sig','gam']
1244        else:
1245            names = ['pos','int','alp','bet','sig','gam']           
1246        head = 13*' '
1247        for name in names:
1248            head += name.center(10)+'esd'.center(10)
1249        print head
1250        if 'C' in dataType:
1251            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1252        else:
1253            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1254        for i,peak in enumerate(Peaks):
1255            ptstr =  ':'
1256            for j in range(len(names)):
1257                name = names[j]
1258                parName = name+str(i)
1259                ptstr += ptfmt[name] % (parmDict[parName])
1260                if parName in varyList:
1261#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1262                    ptstr += ptfmt[name] % (sigDict[parName])
1263                else:
1264#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1265                    ptstr += 10*' '
1266            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1267               
1268    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1269        parmdict.update(zip(varylist,values))
1270        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1271           
1272    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1273        parmdict.update(zip(varylist,values))
1274        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1275        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1276        if dlg:
1277            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1278            if not GoOn:
1279                return -M           #abort!!
1280        return M
1281       
1282    if controls:
1283        Ftol = controls['min dM/M']
1284        derivType = controls['deriv type']
1285    else:
1286        Ftol = 0.0001
1287        derivType = 'analytic'
1288    if oneCycle:
1289        Ftol = 1.0
1290    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1291    yc *= 0.                            #set calcd ones to zero
1292    yb *= 0.
1293    yd *= 0.
1294    xBeg = np.searchsorted(x,Limits[0])
1295    xFin = np.searchsorted(x,Limits[1])
1296    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1297    dataType,insDict,insVary = SetInstParms(Inst)
1298    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1299    parmDict = {}
1300    parmDict.update(bakDict)
1301    parmDict.update(insDict)
1302    parmDict.update(peakDict)
1303    parmDict['Pdabc'] = []      #dummy Pdabc
1304    parmDict.update(Inst2)      #put in real one if there
1305    varyList = bakVary+insVary+peakVary
1306    while True:
1307        begin = time.time()
1308        values =  np.array(Dict2Values(parmDict, varyList))
1309        if FitPgm == 'LSQ':
1310            try:
1311                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1312                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1313                ncyc = int(result[2]['nfev']/2)
1314            finally:
1315                dlg.Destroy()
1316            runtime = time.time()-begin   
1317            chisq = np.sum(result[2]['fvec']**2)
1318            Values2Dict(parmDict, varyList, result[0])
1319            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1320            GOF = chisq/(xFin-xBeg-len(varyList))
1321            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1322            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1323            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1324            try:
1325                sig = np.sqrt(np.diag(result[1])*GOF)
1326                if np.any(np.isnan(sig)):
1327                    print '*** Least squares aborted - some invalid esds possible ***'
1328                break                   #refinement succeeded - finish up!
1329            except ValueError:          #result[1] is None on singular matrix
1330                print '**** Refinement failed - singular matrix ****'
1331                Ipvt = result[2]['ipvt']
1332                for i,ipvt in enumerate(Ipvt):
1333                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1334                        print 'Removing parameter: ',varyList[ipvt-1]
1335                        del(varyList[ipvt-1])
1336                        break
1337        elif FitPgm == 'BFGS':
1338            print 'Other program here'
1339            return
1340       
1341    sigDict = dict(zip(varyList,sig))
1342    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1343    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1344    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1345    GetBackgroundParms(parmDict,Background)
1346    BackgroundPrint(Background,sigDict)
1347    GetInstParms(parmDict,Inst,varyList,Peaks)
1348    InstPrint(Inst,sigDict)
1349    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1350    PeaksPrint(dataType,parmDict,sigDict,varyList)
1351
1352def calcIncident(Iparm,xdata):
1353
1354    def IfunAdv(Iparm,xdata):
1355        Itype = Iparm['Itype']
1356        Icoef = Iparm['Icoeff']
1357        DYI = np.ones((12,xdata.shape[0]))
1358        YI = np.ones_like(xdata)*Icoef[0]
1359       
1360        x = xdata/1000.                 #expressions are in ms
1361        if Itype == 'Exponential':
1362            for i in range(1,10,2):
1363                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1364                YI += Icoef[i]*Eterm
1365                DYI[i] *= Eterm
1366                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1367        elif 'Maxwell'in Itype:
1368            Eterm = np.exp(-Icoef[2]/x**2)
1369            DYI[1] = Eterm/x**5
1370            DYI[2] = -Icoef[1]*DYI[1]/x**2
1371            YI += (Icoef[1]*Eterm/x**5)
1372            if 'Exponential' in Itype:
1373                for i in range(3,12,2):
1374                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1375                    YI += Icoef[i]*Eterm
1376                    DYI[i] *= Eterm
1377                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1378            else:   #Chebyschev
1379                T = (2./x)-1.
1380                Ccof = np.ones((12,xdata.shape[0]))
1381                Ccof[1] = T
1382                for i in range(2,12):
1383                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1384                for i in range(1,10):
1385                    YI += Ccof[i]*Icoef[i+2]
1386                    DYI[i+2] =Ccof[i]
1387        return YI,DYI
1388       
1389    Iesd = np.array(Iparm['Iesd'])
1390    Icovar = Iparm['Icovar']
1391    YI,DYI = IfunAdv(Iparm,xdata)
1392    YI = np.where(YI>0,YI,1.)
1393    WYI = np.zeros_like(xdata)
1394    vcov = np.zeros((12,12))
1395    k = 0
1396    for i in range(12):
1397        for j in range(i,12):
1398            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1399            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1400            k += 1
1401    M = np.inner(vcov,DYI.T)
1402    WYI = np.sum(M*DYI,axis=0)
1403    WYI = np.where(WYI>0.,WYI,0.)
1404    return YI,WYI
1405   
1406#testing data
1407NeedTestData = True
1408def TestData():
1409#    global NeedTestData
1410    NeedTestData = False
1411    global bakType
1412    bakType = 'chebyschev'
1413    global xdata
1414    xdata = np.linspace(4.0,40.0,36000)
1415    global parmDict0
1416    parmDict0 = {
1417        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1418        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1419        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1420        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1421        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1422        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1423        }
1424    global parmDict1
1425    parmDict1 = {
1426        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1427        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1428        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1429        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1430        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1431        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1432        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1433        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1434        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1435        }
1436    global parmDict2
1437    parmDict2 = {
1438        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1439        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1440        'Back0':5.,'Back1':-0.02,'Back2':.004,
1441#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1442        }
1443    global varyList
1444    varyList = []
1445
1446def test0():
1447    if NeedTestData: TestData()
1448    msg = 'test '
1449    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1450    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1451    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1452    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1453    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1454    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1455   
1456def test1():
1457    if NeedTestData: TestData()
1458    time0 = time.time()
1459    for i in range(100):
1460        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1461    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1462   
1463def test2(name,delt):
1464    if NeedTestData: TestData()
1465    varyList = [name,]
1466    xdata = np.linspace(5.6,5.8,400)
1467    hplot = plotter.add('derivatives test for '+name).gca()
1468    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1469    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1470    parmDict2[name] += delt
1471    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1472    hplot.plot(xdata,(y1-y0)/delt,'r+')
1473   
1474def test3(name,delt):
1475    if NeedTestData: TestData()
1476    names = ['pos','sig','gam','shl']
1477    idx = names.index(name)
1478    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1479    xdata = np.linspace(5.6,5.8,800)
1480    dx = xdata[1]-xdata[0]
1481    hplot = plotter.add('derivatives test for '+name).gca()
1482    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1483    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1484    myDict[name] += delt
1485    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1486    hplot.plot(xdata,(y1-y0)/delt,'r+')
1487
1488if __name__ == '__main__':
1489    import GSASIItestplot as plot
1490    global plotter
1491    plotter = plot.PlotNotebook()
1492#    test0()
1493#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1494    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1495        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1496        test2(name,shft)
1497    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1498        test3(name,shft)
1499    print "OK"
1500    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.