source: trunk/GSASIIpwd.py @ 801

Last change on this file since 801 was 801, checked in by vondreele, 9 years ago

fix TOF peak fitting - needed cw*derivs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.2 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2012-11-15 16:44:32 +0000 (Thu, 15 Nov 2012) $
6# $Author: vondreele $
7# $Revision: 801 $
8# $URL: trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 801 2012-11-15 16:44:32Z 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: 801 $")
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                cw = np.diff(xdata[iBeg:iFin])
973                cw = np.append(cw,cw[-1])
974                for i in range(1,6):
975                    dMdpk[i][iBeg:iFin] += intens*cw*dMdipk[i]
976                dMdpk[0][iBeg:iFin] += cw*dMdipk[0]
977                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
978                for parmName in ['pos','int','alp','bet','sig','gam']:
979                    try:
980                        idx = varyList.index(parmName+str(iPeak))
981                        dMdv[idx] = dervDict[parmName]
982                    except ValueError:
983                        pass
984                if 'alpha' in varyList:
985                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
986                if 'beta-0' in varyList:
987                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
988                if 'beta-1' in varyList:
989                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
990                if 'sig-0' in varyList:
991                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
992                if 'sig-1' in varyList:
993                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
994                if 'X' in varyList:
995                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
996                if 'Y' in varyList:
997                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
998                iPeak += 1
999            except KeyError:        #no more peaks to process
1000                break
1001    return dMdv
1002       
1003def Dict2Values(parmdict, varylist):
1004    '''Use before call to leastsq to setup list of values for the parameters
1005    in parmdict, as selected by key in varylist'''
1006    return [parmdict[key] for key in varylist] 
1007   
1008def Values2Dict(parmdict, varylist, values):
1009    ''' Use after call to leastsq to update the parameter dictionary with
1010    values corresponding to keys in varylist'''
1011    parmdict.update(zip(varylist,values))
1012   
1013def SetBackgroundParms(Background):
1014    if len(Background) == 1:            # fix up old backgrounds
1015        BackGround.append({'nDebye':0,'debyeTerms':[]})
1016    bakType,bakFlag = Background[0][:2]
1017    backVals = Background[0][3:]
1018    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1019    Debye = Background[1]           #also has background peaks stuff
1020    backDict = dict(zip(backNames,backVals))
1021    backVary = []
1022    if bakFlag:
1023        backVary = backNames
1024
1025    backDict['nDebye'] = Debye['nDebye']
1026    debyeDict = {}
1027    debyeList = []
1028    for i in range(Debye['nDebye']):
1029        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1030        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1031        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1032    debyeVary = []
1033    for item in debyeList:
1034        if item[1]:
1035            debyeVary.append(item[0])
1036    backDict.update(debyeDict)
1037    backVary += debyeVary
1038
1039    backDict['nPeaks'] = Debye['nPeaks']
1040    peaksDict = {}
1041    peaksList = []
1042    for i in range(Debye['nPeaks']):
1043        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1044        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1045        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1046    peaksVary = []
1047    for item in peaksList:
1048        if item[1]:
1049            peaksVary.append(item[0])
1050    backDict.update(peaksDict)
1051    backVary += peaksVary   
1052    return bakType,backDict,backVary
1053           
1054def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1055   
1056       
1057    def GetBackgroundParms(parmList,Background):
1058        iBak = 0
1059        while True:
1060            try:
1061                bakName = 'Back:'+str(iBak)
1062                Background[0][iBak+3] = parmList[bakName]
1063                iBak += 1
1064            except KeyError:
1065                break
1066        iDb = 0
1067        while True:
1068            names = ['DebyeA:','DebyeR:','DebyeU:']
1069            try:
1070                for i,name in enumerate(names):
1071                    val = parmList[name+str(iDb)]
1072                    Background[1]['debyeTerms'][iDb][2*i] = val
1073                iDb += 1
1074            except KeyError:
1075                break
1076        iDb = 0
1077        while True:
1078            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1079            try:
1080                for i,name in enumerate(names):
1081                    val = parmList[name+str(iDb)]
1082                    Background[1]['peaksList'][iDb][2*i] = val
1083                iDb += 1
1084            except KeyError:
1085                break
1086               
1087    def BackgroundPrint(Background,sigDict):
1088        if Background[0][1]:
1089            print 'Background coefficients for',Background[0][0],'function'
1090            ptfmt = "%12.5f"
1091            ptstr =  'values:'
1092            sigstr = 'esds  :'
1093            for i,back in enumerate(Background[0][3:]):
1094                ptstr += ptfmt % (back)
1095                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1096            print ptstr
1097            print sigstr
1098        else:
1099            print 'Background not refined'
1100        if Background[1]['nDebye']:
1101            parms = ['DebyeA','DebyeR','DebyeU']
1102            print 'Debye diffuse scattering coefficients'
1103            ptfmt = "%12.5f"
1104            names =   'names :'
1105            ptstr =  'values:'
1106            sigstr = 'esds  :'
1107            for item in sigDict:
1108                if 'Debye' in item:
1109                    names += '%12s'%(item)
1110                    sigstr += ptfmt%(sigDict[item])
1111                    parm,id = item.split(':')
1112                    ip = parms.index(parm)
1113                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1114            print names
1115            print ptstr
1116            print sigstr
1117        if Background[1]['nPeaks']:
1118            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1119            print 'Peaks in background coefficients'
1120            ptfmt = "%12.5f"
1121            names =   'names :'
1122            ptstr =  'values:'
1123            sigstr = 'esds  :'
1124            for item in sigDict:
1125                if 'BkPk' in item:
1126                    names += '%12s'%(item)
1127                    sigstr += ptfmt%(sigDict[item])
1128                    parm,id = item.split(':')
1129                    ip = parms.index(parm)
1130                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1131            print names
1132            print ptstr
1133            print sigstr
1134                           
1135    def SetInstParms(Inst):
1136        dataType = Inst['Type'][0]
1137        insVary = []
1138        insNames = []
1139        insVals = []
1140        for parm in Inst:
1141            insNames.append(parm)
1142            insVals.append(Inst[parm][1])
1143            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1144                'beta-0','beta-1','sig-0','sig-1',] and Inst[parm][2]:
1145                    insVary.append(parm)
1146        instDict = dict(zip(insNames,insVals))
1147        instDict['X'] = max(instDict['X'],0.01)
1148        instDict['Y'] = max(instDict['Y'],0.01)
1149        if 'SH/L' in instDict:
1150            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1151        return dataType,instDict,insVary
1152       
1153    def GetInstParms(parmDict,Inst,varyList,Peaks):
1154        for name in Inst:
1155            Inst[name][1] = parmDict[name]
1156        iPeak = 0
1157        while True:
1158            try:
1159                sigName = 'sig'+str(iPeak)
1160                pos = parmDict['pos'+str(iPeak)]
1161                if sigName not in varyList:
1162                    if 'C' in Inst['Type'][0]:
1163                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1164                    else:
1165                        dsp = pos/Inst['difC'][1]
1166                        parmDict[sigName] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1167                gamName = 'gam'+str(iPeak)
1168                if gamName not in varyList:
1169                    if 'C' in Inst['Type'][0]:
1170                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1171                    else:
1172                        dsp = pos/Inst['difC'][1]
1173                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1174                iPeak += 1
1175            except KeyError:
1176                break
1177       
1178    def InstPrint(Inst,sigDict):
1179        print 'Instrument Parameters:'
1180        ptfmt = "%12.6f"
1181        ptlbls = 'names :'
1182        ptstr =  'values:'
1183        sigstr = 'esds  :'
1184        for parm in Inst:
1185            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1186                'beta-0','beta-1','sig-0','sig-1',]:
1187                ptlbls += "%s" % (parm.center(12))
1188                ptstr += ptfmt % (Inst[parm][1])
1189                if parm in sigDict:
1190                    sigstr += ptfmt % (sigDict[parm])
1191                else:
1192                    sigstr += 12*' '
1193        print ptlbls
1194        print ptstr
1195        print sigstr
1196
1197    def SetPeaksParms(dataType,Peaks):
1198        peakNames = []
1199        peakVary = []
1200        peakVals = []
1201        if 'C' in dataType:
1202            names = ['pos','int','sig','gam']
1203        else:
1204            names = ['pos','int','alp','bet','sig','gam']
1205        for i,peak in enumerate(Peaks):
1206            for j,name in enumerate(names):
1207                peakVals.append(peak[2*j])
1208                parName = name+str(i)
1209                peakNames.append(parName)
1210                if peak[2*j+1]:
1211                    peakVary.append(parName)
1212        return dict(zip(peakNames,peakVals)),peakVary
1213               
1214    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1215        if 'C' in Inst['Type'][0]:
1216            names = ['pos','int','sig','gam']
1217        else:
1218            names = ['pos','int','alp','bet','sig','gam']
1219        for i,peak in enumerate(Peaks):
1220            pos = parmDict['pos'+str(i)]
1221            if 'difC' in Inst:
1222                dsp = pos/Inst['difC'][1]
1223            for j in range(len(names)):
1224                parName = names[j]+str(i)
1225                if parName in varyList:
1226                    peak[2*j] = parmDict[parName]
1227                elif 'alpha' in parName:
1228                    peak[2*j] = parmDict['alpha']/dsp
1229                elif 'beta' in parName:
1230                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
1231                elif 'sig' in parName:
1232                    if 'C' in Inst['Type'][0]:
1233                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1234                    else:
1235                        peak[2*j] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1236                elif 'gam' in parName:
1237                    if 'C' in Inst['Type'][0]:
1238                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1239                    else:
1240                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1241                       
1242    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1243        print 'Peak coefficients:'
1244        if 'C' in dataType:
1245            names = ['pos','int','sig','gam']
1246        else:
1247            names = ['pos','int','alp','bet','sig','gam']           
1248        head = 13*' '
1249        for name in names:
1250            head += name.center(10)+'esd'.center(10)
1251        print head
1252        if 'C' in dataType:
1253            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1254        else:
1255            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1256        for i,peak in enumerate(Peaks):
1257            ptstr =  ':'
1258            for j in range(len(names)):
1259                name = names[j]
1260                parName = name+str(i)
1261                ptstr += ptfmt[name] % (parmDict[parName])
1262                if parName in varyList:
1263#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1264                    ptstr += ptfmt[name] % (sigDict[parName])
1265                else:
1266#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1267                    ptstr += 10*' '
1268            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1269               
1270    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1271        parmdict.update(zip(varylist,values))
1272        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1273           
1274    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1275        parmdict.update(zip(varylist,values))
1276        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1277        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1278        if dlg:
1279            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1280            if not GoOn:
1281                return -M           #abort!!
1282        return M
1283       
1284    if controls:
1285        Ftol = controls['min dM/M']
1286        derivType = controls['deriv type']
1287    else:
1288        Ftol = 0.0001
1289        derivType = 'analytic'
1290    if oneCycle:
1291        Ftol = 1.0
1292    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1293    yc *= 0.                            #set calcd ones to zero
1294    yb *= 0.
1295    yd *= 0.
1296    xBeg = np.searchsorted(x,Limits[0])
1297    xFin = np.searchsorted(x,Limits[1])
1298    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1299    dataType,insDict,insVary = SetInstParms(Inst)
1300    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1301    parmDict = {}
1302    parmDict.update(bakDict)
1303    parmDict.update(insDict)
1304    parmDict.update(peakDict)
1305    parmDict['Pdabc'] = []      #dummy Pdabc
1306    parmDict.update(Inst2)      #put in real one if there
1307    varyList = bakVary+insVary+peakVary
1308    while True:
1309        begin = time.time()
1310        values =  np.array(Dict2Values(parmDict, varyList))
1311        if FitPgm == 'LSQ':
1312            try:
1313                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1314                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1315                ncyc = int(result[2]['nfev']/2)
1316            finally:
1317                dlg.Destroy()
1318            runtime = time.time()-begin   
1319            chisq = np.sum(result[2]['fvec']**2)
1320            Values2Dict(parmDict, varyList, result[0])
1321            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1322            GOF = chisq/(xFin-xBeg-len(varyList))
1323            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1324            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1325            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1326            try:
1327                sig = np.sqrt(np.diag(result[1])*GOF)
1328                if np.any(np.isnan(sig)):
1329                    print '*** Least squares aborted - some invalid esds possible ***'
1330                break                   #refinement succeeded - finish up!
1331            except ValueError:          #result[1] is None on singular matrix
1332                print '**** Refinement failed - singular matrix ****'
1333                Ipvt = result[2]['ipvt']
1334                for i,ipvt in enumerate(Ipvt):
1335                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1336                        print 'Removing parameter: ',varyList[ipvt-1]
1337                        del(varyList[ipvt-1])
1338                        break
1339        elif FitPgm == 'BFGS':
1340            print 'Other program here'
1341            return
1342       
1343    sigDict = dict(zip(varyList,sig))
1344    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1345    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1346    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1347    GetBackgroundParms(parmDict,Background)
1348    BackgroundPrint(Background,sigDict)
1349    GetInstParms(parmDict,Inst,varyList,Peaks)
1350    InstPrint(Inst,sigDict)
1351    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1352    PeaksPrint(dataType,parmDict,sigDict,varyList)
1353
1354def calcIncident(Iparm,xdata):
1355
1356    def IfunAdv(Iparm,xdata):
1357        Itype = Iparm['Itype']
1358        Icoef = Iparm['Icoeff']
1359        DYI = np.ones((12,xdata.shape[0]))
1360        YI = np.ones_like(xdata)*Icoef[0]
1361       
1362        x = xdata/1000.                 #expressions are in ms
1363        if Itype == 'Exponential':
1364            for i in range(1,10,2):
1365                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1366                YI += Icoef[i]*Eterm
1367                DYI[i] *= Eterm
1368                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1369        elif 'Maxwell'in Itype:
1370            Eterm = np.exp(-Icoef[2]/x**2)
1371            DYI[1] = Eterm/x**5
1372            DYI[2] = -Icoef[1]*DYI[1]/x**2
1373            YI += (Icoef[1]*Eterm/x**5)
1374            if 'Exponential' in Itype:
1375                for i in range(3,12,2):
1376                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1377                    YI += Icoef[i]*Eterm
1378                    DYI[i] *= Eterm
1379                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1380            else:   #Chebyschev
1381                T = (2./x)-1.
1382                Ccof = np.ones((12,xdata.shape[0]))
1383                Ccof[1] = T
1384                for i in range(2,12):
1385                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1386                for i in range(1,10):
1387                    YI += Ccof[i]*Icoef[i+2]
1388                    DYI[i+2] =Ccof[i]
1389        return YI,DYI
1390       
1391    Iesd = np.array(Iparm['Iesd'])
1392    Icovar = Iparm['Icovar']
1393    YI,DYI = IfunAdv(Iparm,xdata)
1394    YI = np.where(YI>0,YI,1.)
1395    WYI = np.zeros_like(xdata)
1396    vcov = np.zeros((12,12))
1397    k = 0
1398    for i in range(12):
1399        for j in range(i,12):
1400            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1401            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1402            k += 1
1403    M = np.inner(vcov,DYI.T)
1404    WYI = np.sum(M*DYI,axis=0)
1405    WYI = np.where(WYI>0.,WYI,0.)
1406    return YI,WYI
1407   
1408#testing data
1409NeedTestData = True
1410def TestData():
1411#    global NeedTestData
1412    NeedTestData = False
1413    global bakType
1414    bakType = 'chebyschev'
1415    global xdata
1416    xdata = np.linspace(4.0,40.0,36000)
1417    global parmDict0
1418    parmDict0 = {
1419        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1420        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1421        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1422        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1423        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1424        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1425        }
1426    global parmDict1
1427    parmDict1 = {
1428        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1429        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1430        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1431        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1432        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1433        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1434        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1435        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1436        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1437        }
1438    global parmDict2
1439    parmDict2 = {
1440        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1441        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1442        'Back0':5.,'Back1':-0.02,'Back2':.004,
1443#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1444        }
1445    global varyList
1446    varyList = []
1447
1448def test0():
1449    if NeedTestData: TestData()
1450    msg = 'test '
1451    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1452    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1453    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1454    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1455    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1456    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1457   
1458def test1():
1459    if NeedTestData: TestData()
1460    time0 = time.time()
1461    for i in range(100):
1462        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1463    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1464   
1465def test2(name,delt):
1466    if NeedTestData: TestData()
1467    varyList = [name,]
1468    xdata = np.linspace(5.6,5.8,400)
1469    hplot = plotter.add('derivatives test for '+name).gca()
1470    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1471    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1472    parmDict2[name] += delt
1473    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1474    hplot.plot(xdata,(y1-y0)/delt,'r+')
1475   
1476def test3(name,delt):
1477    if NeedTestData: TestData()
1478    names = ['pos','sig','gam','shl']
1479    idx = names.index(name)
1480    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1481    xdata = np.linspace(5.6,5.8,800)
1482    dx = xdata[1]-xdata[0]
1483    hplot = plotter.add('derivatives test for '+name).gca()
1484    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1485    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1486    myDict[name] += delt
1487    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1488    hplot.plot(xdata,(y1-y0)/delt,'r+')
1489
1490if __name__ == '__main__':
1491    import GSASIItestplot as plot
1492    global plotter
1493    plotter = plot.PlotNotebook()
1494#    test0()
1495#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1496    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1497        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1498        test2(name,shft)
1499    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1500        test3(name,shft)
1501    print "OK"
1502    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.