source: trunk/GSASIIpwd.py @ 803

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

GSASIImath.py: map search gives mag in name, & change a looping form in findOffset
GSASIIphsGUI.py: map search gives mag in name
GSASIIplot.py: more curves for instrument parms plot
GSASIIpwd.py: make cw array of channel widths & use them
GSASIIstruct.py: ditto
G2pwd_fxye.py: fix for zero sig

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 56.8 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3#GSASII powder calculation module
4########### SVN repository information ###################
5# $Date: 2012-11-21 19:30:38 +0000 (Wed, 21 Nov 2012) $
6# $Author: vondreele $
7# $Revision: 803 $
8# $URL: trunk/GSASIIpwd.py $
9# $Id: GSASIIpwd.py 803 2012-11-21 19:30:38Z 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: 803 $")
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    cw = np.diff(xdata)
453    cw = np.append(cw,cw[-1])
454    while True:
455        key = pfx+'Back:'+str(nBak)
456        if key in parmDict:
457            nBak += 1
458        else:
459            break
460    if bakType in ['chebyschev','cosine']:
461        for iBak in range(nBak):   
462            key = pfx+'Back:'+str(iBak)
463            if bakType == 'chebyschev':
464                yb += parmDict[key]*(xdata-xdata[0])**iBak
465            elif bakType == 'cosine':
466                yb += parmDict[key]*npcosd(xdata*iBak)
467    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
468        if nBak == 1:
469            yb = np.ones_like(xdata)*parmDict[pfx+'Back:0']
470        elif nBak == 2:
471            dX = xdata[-1]-xdata[0]
472            T2 = (xdata-xdata[0])/dX
473            T1 = 1.0-T2
474            yb = parmDict[pfx+'Back:0']*T1+parmDict[pfx+'Back:1']*T2
475        else:
476            if bakType == 'lin interpolate':
477                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
478            elif bakType == 'inv interpolate':
479                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
480            elif bakType == 'log interpolate':
481                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
482            bakPos[0] = xdata[0]
483            bakPos[-1] = xdata[-1]
484            bakVals = np.zeros(nBak)
485            for i in range(nBak):
486                bakVals[i] = parmDict[pfx+'Back:'+str(i)]
487            bakInt = si.interp1d(bakPos,bakVals,'linear')
488            yb = bakInt(xdata)
489    if 'difC' in parmDict:
490        ff = 1.
491    else:       
492        try:
493            wave = parmDict[pfx+'Lam']
494        except KeyError:
495            wave = parmDict[pfx+'Lam1']
496        q = 4.0*np.pi*npsind(xdata/2.0)/wave
497        SQ = (q/(4.*np.pi))**2
498        FF = G2elem.GetFormFactorCoeff('Si')[0]
499        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
500    iD = 0       
501    while True:
502        try:
503            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
504            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
505            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
506            yb += ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
507            iD += 1       
508        except KeyError:
509            break
510    iD = 0
511    while True:
512        try:
513            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
514            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
515            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
516            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
517            shl = 0.002
518            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
519            iBeg = np.searchsorted(xdata,pkP-fmin)
520            iFin = np.searchsorted(xdata,pkP+fmax)
521            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
522            iD += 1       
523        except KeyError:
524            break       
525    return yb
526   
527def getBackgroundDerv(pfx,parmDict,bakType,xdata):
528    nBak = 0
529    while True:
530        key = pfx+'Back:'+str(nBak)
531        if key in parmDict:
532            nBak += 1
533        else:
534            break
535    dydb = np.zeros(shape=(nBak,len(xdata)))
536    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
537    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
538    cw = np.diff(xdata)
539    cw = np.append(cw,cw[-1])
540
541    if bakType in ['chebyschev','cosine']:
542        for iBak in range(nBak):   
543            if bakType == 'chebyschev':
544                dydb[iBak] = (xdata-xdata[0])**iBak
545            elif bakType == 'cosine':
546                dydb[iBak] = npcosd(xdata*iBak)
547    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
548        if nBak == 1:
549            dydb[0] = np.ones_like(xdata)
550        elif nBak == 2:
551            dX = xdata[-1]-xdata[0]
552            T2 = (xdata-xdata[0])/dX
553            T1 = 1.0-T2
554            dydb = [T1,T2]
555        else:
556            if bakType == 'lin interpolate':
557                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
558            elif bakType == 'inv interpolate':
559                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
560            elif bakType == 'log interpolate':
561                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
562            bakPos[0] = xdata[0]
563            bakPos[-1] = xdata[-1]
564            for i,pos in enumerate(bakPos):
565                if i == 0:
566                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
567                elif i == len(bakPos)-1:
568                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
569                else:
570                    dydb[i] = np.where(xdata>bakPos[i],
571                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
572                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
573    if 'difC' in parmDict:
574        ff = 1.
575    else:
576        try:
577            wave = parmDict[pfx+'Lam']
578        except KeyError:
579            wave = parmDict[pfx+'Lam1']
580        q = 4.0*np.pi*npsind(xdata/2.0)/wave
581        SQ = (q/(4*np.pi))**2
582        FF = G2elem.GetFormFactorCoeff('Si')[0]
583        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
584    iD = 0       
585    while True:
586        try:
587            dbA = parmDict[pfx+'DebyeA:'+str(iD)]
588            dbR = parmDict[pfx+'DebyeR:'+str(iD)]
589            dbU = parmDict[pfx+'DebyeU:'+str(iD)]
590            sqr = np.sin(q*dbR)/(q*dbR)
591            cqr = np.cos(q*dbR)
592            temp = np.exp(-dbU*q**2)
593            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
594            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
595            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
596            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
597        except KeyError:
598            break
599    iD = 0
600    while True:
601        try:
602            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
603            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
604            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
605            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
606            shl = 0.002
607            Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,shl)
608            iBeg = np.searchsorted(xdata,pkP-fmin)
609            iFin = np.searchsorted(xdata,pkP+fmax)
610            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
611            dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
612            dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
613            dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
614            dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
615            iD += 1       
616        except KeyError:
617            break
618    return dydb,dyddb,dydpk
619
620#use old fortran routine
621def getFCJVoigt3(pos,sig,gam,shl,xdata):
622   
623    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
624#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
625    Df /= np.sum(Df)
626    return Df
627
628def getdFCJVoigt3(pos,sig,gam,shl,xdata):
629   
630    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
631#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
632    sumDf = np.sum(Df)
633    return Df,dFdp,dFds,dFdg,dFdsh
634
635def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
636    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
637    Df /= np.sum(Df)
638    return Df 
639   
640def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
641    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
642    sumDf = np.sum(Df)
643    return Df,dFdp,dFda,dFdb,dFds,dFdg   
644
645def ellipseSize(H,Sij,GB):
646    HX = np.inner(H.T,GB)
647    lenHX = np.sqrt(np.sum(HX**2))
648    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
649    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
650    lenR = np.sqrt(np.sum(R**2))
651    return lenR
652
653def ellipseSizeDerv(H,Sij,GB):
654    lenR = ellipseSize(H,Sij,GB)
655    delt = 0.001
656    dRdS = np.zeros(6)
657    for i in range(6):
658        dSij = Sij[:]
659        dSij[i] += delt
660        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
661    return lenR,dRdS
662
663def getHKLpeak(dmin,SGData,A):
664    HKL = G2lat.GenHLaue(dmin,SGData,A)       
665    HKLs = []
666    for h,k,l,d in HKL:
667        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
668        if not ext:
669            HKLs.append([h,k,l,d,-1])
670    return HKLs
671
672def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
673   
674    yb = getBackground('',parmDict,bakType,xdata)
675    yc = np.zeros_like(yb)
676    cw = np.diff(xdata)
677    cw = np.append(cw,cw[-1])
678    if 'C' in dataType:
679        U = parmDict['U']
680        V = parmDict['V']
681        W = parmDict['W']
682        X = parmDict['X']
683        Y = parmDict['Y']
684        shl = max(parmDict['SH/L'],0.002)
685        Ka2 = False
686        if 'Lam1' in parmDict.keys():
687            Ka2 = True
688            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
689            kRatio = parmDict['I(L2)/I(L1)']
690        iPeak = 0
691        while True:
692            try:
693                pos = parmDict['pos'+str(iPeak)]
694                theta = (pos-parmDict['Zero'])/2.0
695                intens = parmDict['int'+str(iPeak)]
696                sigName = 'sig'+str(iPeak)
697                if sigName in varyList:
698                    sig = parmDict[sigName]
699                else:
700                    sig = U*tand(theta)**2+V*tand(theta)+W
701                sig = max(sig,0.001)          #avoid neg sigma
702                gamName = 'gam'+str(iPeak)
703                if gamName in varyList:
704                    gam = parmDict[gamName]
705                else:
706                    gam = X/cosd(theta)+Y*tand(theta)
707                gam = max(gam,0.001)             #avoid neg gamma
708                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
709                iBeg = np.searchsorted(xdata,pos-fmin)
710                iFin = np.searchsorted(xdata,pos+fmin)
711                if not iBeg+iFin:       #peak below low limit
712                    iPeak += 1
713                    continue
714                elif not iBeg-iFin:     #peak above high limit
715                    return yb+yc
716                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
717                if Ka2:
718                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
719                    iBeg = np.searchsorted(xdata,pos2-fmin)
720                    iFin = np.searchsorted(xdata,pos2+fmin)
721                    if iBeg-iFin:
722                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
723                iPeak += 1
724            except KeyError:        #no more peaks to process
725                return yb+yc
726    else:
727        Pdabc = parmDict['Pdabc']
728        difC = parmDict['difC']
729        alp0 = parmDict['alpha']
730        bet0 = parmDict['beta-0']
731        bet1 = parmDict['beta-1']
732        sig0 = parmDict['sig-0']
733        sig1 = parmDict['sig-1']
734        X = parmDict['X']
735        Y = parmDict['Y']
736        iPeak = 0
737        while True:
738            try:
739                pos = parmDict['pos'+str(iPeak)]               
740                tof = pos-parmDict['Zero']
741                dsp = tof/difC
742                intens = parmDict['int'+str(iPeak)]
743                alpName = 'alp'+str(iPeak)
744                if alpName in varyList:
745                    alp = parmDict[alpName]
746                else:
747                    if len(Pdabc):
748                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
749                    else:
750                        alp = alp0/dsp
751                betName = 'bet'+str(iPeak)
752                if betName in varyList:
753                    bet = parmDict[betName]
754                else:
755                    if len(Pdabc):
756                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
757                    else:
758                        bet = bet0+bet1/dsp**4
759                sigName = 'sig'+str(iPeak)
760                if sigName in varyList:
761                    sig = parmDict[sigName]
762                else:
763                    sig = sig0+sig1*dsp**2
764                gamName = 'gam'+str(iPeak)
765                if gamName in varyList:
766                    gam = parmDict[gamName]
767                else:
768                    gam = X*dsp**2+Y*dsp
769                gam = max(gam,0.001)             #avoid neg gamma
770                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
771                iBeg = np.searchsorted(xdata,pos-fmin)
772                iFin = np.searchsorted(xdata,pos+fmax)
773                lenX = len(xdata)
774                if not iBeg:
775                    iFin = np.searchsorted(xdata,pos+fmax)
776                elif iBeg == lenX:
777                    iFin = iBeg
778                else:
779                    iFin = np.searchsorted(xdata,pos+fmax)
780                if not iBeg+iFin:       #peak below low limit
781                    iPeak += 1
782                    continue
783                elif not iBeg-iFin:     #peak above high limit
784                    return yb+yc
785                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
786                iPeak += 1
787            except KeyError:        #no more peaks to process
788                return yb+yc
789           
790def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
791# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
792    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
793    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
794    if 'Back:0' in varyList:            #background derivs are in front if present
795        dMdv[0:len(dMdb)] = dMdb
796    names = ['DebyeA','DebyeR','DebyeU']
797    for name in varyList:
798        if 'Debye' in name:
799            parm,id = name.split(':')
800            ip = names.index(parm)
801            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
802    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
803    for name in varyList:
804        if 'BkPk' in name:
805            parm,id = name.split(':')
806            ip = names.index(parm)
807            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
808    cw = np.diff(xdata)
809    cw = np.append(cw,cw[-1])
810    if 'C' in dataType:
811        U = parmDict['U']
812        V = parmDict['V']
813        W = parmDict['W']
814        X = parmDict['X']
815        Y = parmDict['Y']
816        shl = max(parmDict['SH/L'],0.002)
817        Ka2 = False
818        if 'Lam1' in parmDict.keys():
819            Ka2 = True
820            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
821            kRatio = parmDict['I(L2)/I(L1)']
822        iPeak = 0
823        while True:
824            try:
825                pos = parmDict['pos'+str(iPeak)]
826                theta = (pos-parmDict['Zero'])/2.0
827                intens = parmDict['int'+str(iPeak)]
828                sigName = 'sig'+str(iPeak)
829                tanth = tand(theta)
830                costh = cosd(theta)
831                if sigName in varyList:
832                    sig = parmDict[sigName]
833                else:
834                    sig = U*tanth**2+V*tanth+W
835                    dsdU = tanth**2
836                    dsdV = tanth
837                    dsdW = 1.0
838                sig = max(sig,0.001)          #avoid neg sigma
839                gamName = 'gam'+str(iPeak)
840                if gamName in varyList:
841                    gam = parmDict[gamName]
842                else:
843                    gam = X/costh+Y*tanth
844                    dgdX = 1.0/costh
845                    dgdY = tanth
846                gam = max(gam,0.001)             #avoid neg gamma
847                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
848                iBeg = np.searchsorted(xdata,pos-fmin)
849                iFin = np.searchsorted(xdata,pos+fmin)
850                if not iBeg+iFin:       #peak below low limit
851                    iPeak += 1
852                    continue
853                elif not iBeg-iFin:     #peak above high limit
854                    break
855                dMdpk = np.zeros(shape=(6,len(xdata)))
856                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
857                for i in range(1,5):
858                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
859                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
860                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
861                if Ka2:
862                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
863                    iBeg = np.searchsorted(xdata,pos2-fmin)
864                    iFin = np.searchsorted(xdata,pos2+fmin)
865                    if iBeg-iFin:
866                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
867                        for i in range(1,5):
868                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
869                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
870                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
871                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
872                for parmName in ['pos','int','sig','gam']:
873                    try:
874                        idx = varyList.index(parmName+str(iPeak))
875                        dMdv[idx] = dervDict[parmName]
876                    except ValueError:
877                        pass
878                if 'U' in varyList:
879                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
880                if 'V' in varyList:
881                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
882                if 'W' in varyList:
883                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
884                if 'X' in varyList:
885                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
886                if 'Y' in varyList:
887                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
888                if 'SH/L' in varyList:
889                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
890                if 'I(L2)/I(L1)' in varyList:
891                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
892                iPeak += 1
893            except KeyError:        #no more peaks to process
894                break
895    else:
896        Pdabc = parmDict['Pdabc']
897        difC = parmDict['difC']
898        alp0 = parmDict['alpha']
899        bet0 = parmDict['beta-0']
900        bet1 = parmDict['beta-1']
901        sig0 = parmDict['sig-0']
902        sig1 = parmDict['sig-1']
903        X = parmDict['X']
904        Y = parmDict['Y']
905        iPeak = 0
906        while True:
907            try:
908                pos = parmDict['pos'+str(iPeak)]               
909                tof = pos-parmDict['Zero']
910                dsp = tof/difC
911                intens = parmDict['int'+str(iPeak)]
912                alpName = 'alp'+str(iPeak)
913                if alpName in varyList:
914                    alp = parmDict[alpName]
915                else:
916                    if len(Pdabc):
917                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
918                    else:
919                        alp = alp0/dsp
920                    dada0 = 1./dsp
921                betName = 'bet'+str(iPeak)
922                if betName in varyList:
923                    bet = parmDict[betName]
924                else:
925                    if len(Pdabc):
926                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
927                    else:
928                        bet = bet0+bet1/dsp**4
929                    dbdb0 = 1.
930                    dbdb1 = 1/dsp**4
931                sigName = 'sig'+str(iPeak)
932                if sigName in varyList:
933                    sig = parmDict[sigName]
934                else:
935                    sig = sig0+sig1*dsp**2
936                    dsds0 = 1.
937                    dsds1 = dsp**2
938                gamName = 'gam'+str(iPeak)
939                if gamName in varyList:
940                    gam = parmDict[gamName]
941                else:
942                    gam = X*dsp**2+Y*dsp
943                    dsdX = dsp**2
944                    dsdY = dsp
945                gam = max(gam,0.001)             #avoid neg gamma
946                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
947                iBeg = np.searchsorted(xdata,pos-fmin)
948                lenX = len(xdata)
949                if not iBeg:
950                    iFin = np.searchsorted(xdata,pos+fmax)
951                elif iBeg == lenX:
952                    iFin = iBeg
953                else:
954                    iFin = np.searchsorted(xdata,pos+fmax)
955                if not iBeg+iFin:       #peak below low limit
956                    iPeak += 1
957                    continue
958                elif not iBeg-iFin:     #peak above high limit
959                    break
960                dMdpk = np.zeros(shape=(7,len(xdata)))
961                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
962                for i in range(1,6):
963                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
964                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
965                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
966                for parmName in ['pos','int','alp','bet','sig','gam']:
967                    try:
968                        idx = varyList.index(parmName+str(iPeak))
969                        dMdv[idx] = dervDict[parmName]
970                    except ValueError:
971                        pass
972                if 'alpha' in varyList:
973                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
974                if 'beta-0' in varyList:
975                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
976                if 'beta-1' in varyList:
977                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
978                if 'sig-0' in varyList:
979                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
980                if 'sig-1' in varyList:
981                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
982                if 'X' in varyList:
983                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
984                if 'Y' in varyList:
985                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
986                iPeak += 1
987            except KeyError:        #no more peaks to process
988                break
989    return dMdv
990       
991def Dict2Values(parmdict, varylist):
992    '''Use before call to leastsq to setup list of values for the parameters
993    in parmdict, as selected by key in varylist'''
994    return [parmdict[key] for key in varylist] 
995   
996def Values2Dict(parmdict, varylist, values):
997    ''' Use after call to leastsq to update the parameter dictionary with
998    values corresponding to keys in varylist'''
999    parmdict.update(zip(varylist,values))
1000   
1001def SetBackgroundParms(Background):
1002    if len(Background) == 1:            # fix up old backgrounds
1003        BackGround.append({'nDebye':0,'debyeTerms':[]})
1004    bakType,bakFlag = Background[0][:2]
1005    backVals = Background[0][3:]
1006    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1007    Debye = Background[1]           #also has background peaks stuff
1008    backDict = dict(zip(backNames,backVals))
1009    backVary = []
1010    if bakFlag:
1011        backVary = backNames
1012
1013    backDict['nDebye'] = Debye['nDebye']
1014    debyeDict = {}
1015    debyeList = []
1016    for i in range(Debye['nDebye']):
1017        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1018        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1019        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1020    debyeVary = []
1021    for item in debyeList:
1022        if item[1]:
1023            debyeVary.append(item[0])
1024    backDict.update(debyeDict)
1025    backVary += debyeVary
1026
1027    backDict['nPeaks'] = Debye['nPeaks']
1028    peaksDict = {}
1029    peaksList = []
1030    for i in range(Debye['nPeaks']):
1031        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1032        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1033        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1034    peaksVary = []
1035    for item in peaksList:
1036        if item[1]:
1037            peaksVary.append(item[0])
1038    backDict.update(peaksDict)
1039    backVary += peaksVary   
1040    return bakType,backDict,backVary
1041           
1042def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1043   
1044       
1045    def GetBackgroundParms(parmList,Background):
1046        iBak = 0
1047        while True:
1048            try:
1049                bakName = 'Back:'+str(iBak)
1050                Background[0][iBak+3] = parmList[bakName]
1051                iBak += 1
1052            except KeyError:
1053                break
1054        iDb = 0
1055        while True:
1056            names = ['DebyeA:','DebyeR:','DebyeU:']
1057            try:
1058                for i,name in enumerate(names):
1059                    val = parmList[name+str(iDb)]
1060                    Background[1]['debyeTerms'][iDb][2*i] = val
1061                iDb += 1
1062            except KeyError:
1063                break
1064        iDb = 0
1065        while True:
1066            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1067            try:
1068                for i,name in enumerate(names):
1069                    val = parmList[name+str(iDb)]
1070                    Background[1]['peaksList'][iDb][2*i] = val
1071                iDb += 1
1072            except KeyError:
1073                break
1074               
1075    def BackgroundPrint(Background,sigDict):
1076        if Background[0][1]:
1077            print 'Background coefficients for',Background[0][0],'function'
1078            ptfmt = "%12.5f"
1079            ptstr =  'values:'
1080            sigstr = 'esds  :'
1081            for i,back in enumerate(Background[0][3:]):
1082                ptstr += ptfmt % (back)
1083                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1084            print ptstr
1085            print sigstr
1086        else:
1087            print 'Background not refined'
1088        if Background[1]['nDebye']:
1089            parms = ['DebyeA','DebyeR','DebyeU']
1090            print 'Debye diffuse scattering coefficients'
1091            ptfmt = "%12.5f"
1092            names =   'names :'
1093            ptstr =  'values:'
1094            sigstr = 'esds  :'
1095            for item in sigDict:
1096                if 'Debye' in item:
1097                    names += '%12s'%(item)
1098                    sigstr += ptfmt%(sigDict[item])
1099                    parm,id = item.split(':')
1100                    ip = parms.index(parm)
1101                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1102            print names
1103            print ptstr
1104            print sigstr
1105        if Background[1]['nPeaks']:
1106            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1107            print 'Peaks in background coefficients'
1108            ptfmt = "%12.5f"
1109            names =   'names :'
1110            ptstr =  'values:'
1111            sigstr = 'esds  :'
1112            for item in sigDict:
1113                if 'BkPk' in item:
1114                    names += '%12s'%(item)
1115                    sigstr += ptfmt%(sigDict[item])
1116                    parm,id = item.split(':')
1117                    ip = parms.index(parm)
1118                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1119            print names
1120            print ptstr
1121            print sigstr
1122                           
1123    def SetInstParms(Inst):
1124        dataType = Inst['Type'][0]
1125        insVary = []
1126        insNames = []
1127        insVals = []
1128        for parm in Inst:
1129            insNames.append(parm)
1130            insVals.append(Inst[parm][1])
1131            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1132                'beta-0','beta-1','sig-0','sig-1',] and Inst[parm][2]:
1133                    insVary.append(parm)
1134        instDict = dict(zip(insNames,insVals))
1135        instDict['X'] = max(instDict['X'],0.01)
1136        instDict['Y'] = max(instDict['Y'],0.01)
1137        if 'SH/L' in instDict:
1138            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1139        return dataType,instDict,insVary
1140       
1141    def GetInstParms(parmDict,Inst,varyList,Peaks):
1142        for name in Inst:
1143            Inst[name][1] = parmDict[name]
1144        iPeak = 0
1145        while True:
1146            try:
1147                sigName = 'sig'+str(iPeak)
1148                pos = parmDict['pos'+str(iPeak)]
1149                if sigName not in varyList:
1150                    if 'C' in Inst['Type'][0]:
1151                        parmDict[sigName] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1152                    else:
1153                        dsp = pos/Inst['difC'][1]
1154                        parmDict[sigName] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1155                gamName = 'gam'+str(iPeak)
1156                if gamName not in varyList:
1157                    if 'C' in Inst['Type'][0]:
1158                        parmDict[gamName] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1159                    else:
1160                        dsp = pos/Inst['difC'][1]
1161                        parmDict[sigName] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1162                iPeak += 1
1163            except KeyError:
1164                break
1165       
1166    def InstPrint(Inst,sigDict):
1167        print 'Instrument Parameters:'
1168        ptfmt = "%12.6f"
1169        ptlbls = 'names :'
1170        ptstr =  'values:'
1171        sigstr = 'esds  :'
1172        for parm in Inst:
1173            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1174                'beta-0','beta-1','sig-0','sig-1',]:
1175                ptlbls += "%s" % (parm.center(12))
1176                ptstr += ptfmt % (Inst[parm][1])
1177                if parm in sigDict:
1178                    sigstr += ptfmt % (sigDict[parm])
1179                else:
1180                    sigstr += 12*' '
1181        print ptlbls
1182        print ptstr
1183        print sigstr
1184
1185    def SetPeaksParms(dataType,Peaks):
1186        peakNames = []
1187        peakVary = []
1188        peakVals = []
1189        if 'C' in dataType:
1190            names = ['pos','int','sig','gam']
1191        else:
1192            names = ['pos','int','alp','bet','sig','gam']
1193        for i,peak in enumerate(Peaks):
1194            for j,name in enumerate(names):
1195                peakVals.append(peak[2*j])
1196                parName = name+str(i)
1197                peakNames.append(parName)
1198                if peak[2*j+1]:
1199                    peakVary.append(parName)
1200        return dict(zip(peakNames,peakVals)),peakVary
1201               
1202    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1203        if 'C' in Inst['Type'][0]:
1204            names = ['pos','int','sig','gam']
1205        else:
1206            names = ['pos','int','alp','bet','sig','gam']
1207        for i,peak in enumerate(Peaks):
1208            pos = parmDict['pos'+str(i)]
1209            if 'difC' in Inst:
1210                dsp = pos/Inst['difC'][1]
1211            for j in range(len(names)):
1212                parName = names[j]+str(i)
1213                if parName in varyList:
1214                    peak[2*j] = parmDict[parName]
1215                elif 'alpha' in parName:
1216                    peak[2*j] = parmDict['alpha']/dsp
1217                elif 'beta' in parName:
1218                    peak[2*j] = parmDict['beta-0']+parmDict['beta-1']/dsp**4
1219                elif 'sig' in parName:
1220                    if 'C' in Inst['Type'][0]:
1221                        peak[2*j] = parmDict['U']*tand(pos/2.0)**2+parmDict['V']*tand(pos/2.0)+parmDict['W']
1222                    else:
1223                        peak[2*j] = parmDict['sig-0']+parmDict['sig-1']*dsp**2
1224                elif 'gam' in parName:
1225                    if 'C' in Inst['Type'][0]:
1226                        peak[2*j] = parmDict['X']/cosd(pos/2.0)+parmDict['Y']*tand(pos/2.0)
1227                    else:
1228                        peak[2*j] = parmDict['X']*dsp**2+parmDict['Y']*dsp
1229                       
1230    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1231        print 'Peak coefficients:'
1232        if 'C' in dataType:
1233            names = ['pos','int','sig','gam']
1234        else:
1235            names = ['pos','int','alp','bet','sig','gam']           
1236        head = 13*' '
1237        for name in names:
1238            head += name.center(10)+'esd'.center(10)
1239        print head
1240        if 'C' in dataType:
1241            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1242        else:
1243            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1244        for i,peak in enumerate(Peaks):
1245            ptstr =  ':'
1246            for j in range(len(names)):
1247                name = names[j]
1248                parName = name+str(i)
1249                ptstr += ptfmt[name] % (parmDict[parName])
1250                if parName in varyList:
1251#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1252                    ptstr += ptfmt[name] % (sigDict[parName])
1253                else:
1254#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1255                    ptstr += 10*' '
1256            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1257               
1258    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1259        parmdict.update(zip(varylist,values))
1260        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1261           
1262    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1263        parmdict.update(zip(varylist,values))
1264        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1265        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1266        if dlg:
1267            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1268            if not GoOn:
1269                return -M           #abort!!
1270        return M
1271       
1272    if controls:
1273        Ftol = controls['min dM/M']
1274        derivType = controls['deriv type']
1275    else:
1276        Ftol = 0.0001
1277        derivType = 'analytic'
1278    if oneCycle:
1279        Ftol = 1.0
1280    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1281    yc *= 0.                            #set calcd ones to zero
1282    yb *= 0.
1283    yd *= 0.
1284    xBeg = np.searchsorted(x,Limits[0])
1285    xFin = np.searchsorted(x,Limits[1])
1286    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1287    dataType,insDict,insVary = SetInstParms(Inst)
1288    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1289    parmDict = {}
1290    parmDict.update(bakDict)
1291    parmDict.update(insDict)
1292    parmDict.update(peakDict)
1293    parmDict['Pdabc'] = []      #dummy Pdabc
1294    parmDict.update(Inst2)      #put in real one if there
1295    varyList = bakVary+insVary+peakVary
1296    while True:
1297        begin = time.time()
1298        values =  np.array(Dict2Values(parmDict, varyList))
1299        if FitPgm == 'LSQ':
1300            try:
1301                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1302                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1303                ncyc = int(result[2]['nfev']/2)
1304            finally:
1305                dlg.Destroy()
1306            runtime = time.time()-begin   
1307            chisq = np.sum(result[2]['fvec']**2)
1308            Values2Dict(parmDict, varyList, result[0])
1309            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1310            GOF = chisq/(xFin-xBeg-len(varyList))
1311            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1312            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1313            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1314            try:
1315                sig = np.sqrt(np.diag(result[1])*GOF)
1316                if np.any(np.isnan(sig)):
1317                    print '*** Least squares aborted - some invalid esds possible ***'
1318                break                   #refinement succeeded - finish up!
1319            except ValueError:          #result[1] is None on singular matrix
1320                print '**** Refinement failed - singular matrix ****'
1321                Ipvt = result[2]['ipvt']
1322                for i,ipvt in enumerate(Ipvt):
1323                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1324                        print 'Removing parameter: ',varyList[ipvt-1]
1325                        del(varyList[ipvt-1])
1326                        break
1327        elif FitPgm == 'BFGS':
1328            print 'Other program here'
1329            return
1330       
1331    sigDict = dict(zip(varyList,sig))
1332    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1333    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1334    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1335    GetBackgroundParms(parmDict,Background)
1336    BackgroundPrint(Background,sigDict)
1337    GetInstParms(parmDict,Inst,varyList,Peaks)
1338    InstPrint(Inst,sigDict)
1339    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1340    PeaksPrint(dataType,parmDict,sigDict,varyList)
1341
1342def calcIncident(Iparm,xdata):
1343
1344    def IfunAdv(Iparm,xdata):
1345        Itype = Iparm['Itype']
1346        Icoef = Iparm['Icoeff']
1347        DYI = np.ones((12,xdata.shape[0]))
1348        YI = np.ones_like(xdata)*Icoef[0]
1349       
1350        x = xdata/1000.                 #expressions are in ms
1351        if Itype == 'Exponential':
1352            for i in range(1,10,2):
1353                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1354                YI += Icoef[i]*Eterm
1355                DYI[i] *= Eterm
1356                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1357        elif 'Maxwell'in Itype:
1358            Eterm = np.exp(-Icoef[2]/x**2)
1359            DYI[1] = Eterm/x**5
1360            DYI[2] = -Icoef[1]*DYI[1]/x**2
1361            YI += (Icoef[1]*Eterm/x**5)
1362            if 'Exponential' in Itype:
1363                for i in range(3,12,2):
1364                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1365                    YI += Icoef[i]*Eterm
1366                    DYI[i] *= Eterm
1367                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1368            else:   #Chebyschev
1369                T = (2./x)-1.
1370                Ccof = np.ones((12,xdata.shape[0]))
1371                Ccof[1] = T
1372                for i in range(2,12):
1373                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1374                for i in range(1,10):
1375                    YI += Ccof[i]*Icoef[i+2]
1376                    DYI[i+2] =Ccof[i]
1377        return YI,DYI
1378       
1379    Iesd = np.array(Iparm['Iesd'])
1380    Icovar = Iparm['Icovar']
1381    YI,DYI = IfunAdv(Iparm,xdata)
1382    YI = np.where(YI>0,YI,1.)
1383    WYI = np.zeros_like(xdata)
1384    vcov = np.zeros((12,12))
1385    k = 0
1386    for i in range(12):
1387        for j in range(i,12):
1388            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1389            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1390            k += 1
1391    M = np.inner(vcov,DYI.T)
1392    WYI = np.sum(M*DYI,axis=0)
1393    WYI = np.where(WYI>0.,WYI,0.)
1394    return YI,WYI
1395   
1396#testing data
1397NeedTestData = True
1398def TestData():
1399#    global NeedTestData
1400    NeedTestData = False
1401    global bakType
1402    bakType = 'chebyschev'
1403    global xdata
1404    xdata = np.linspace(4.0,40.0,36000)
1405    global parmDict0
1406    parmDict0 = {
1407        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1408        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1409        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1410        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1411        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1412        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1413        }
1414    global parmDict1
1415    parmDict1 = {
1416        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1417        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1418        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1419        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1420        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1421        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1422        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1423        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1424        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1425        }
1426    global parmDict2
1427    parmDict2 = {
1428        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1429        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1430        'Back0':5.,'Back1':-0.02,'Back2':.004,
1431#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1432        }
1433    global varyList
1434    varyList = []
1435
1436def test0():
1437    if NeedTestData: TestData()
1438    msg = 'test '
1439    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1440    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1441    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1442    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1443    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1444    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1445   
1446def test1():
1447    if NeedTestData: TestData()
1448    time0 = time.time()
1449    for i in range(100):
1450        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1451    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1452   
1453def test2(name,delt):
1454    if NeedTestData: TestData()
1455    varyList = [name,]
1456    xdata = np.linspace(5.6,5.8,400)
1457    hplot = plotter.add('derivatives test for '+name).gca()
1458    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1459    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1460    parmDict2[name] += delt
1461    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1462    hplot.plot(xdata,(y1-y0)/delt,'r+')
1463   
1464def test3(name,delt):
1465    if NeedTestData: TestData()
1466    names = ['pos','sig','gam','shl']
1467    idx = names.index(name)
1468    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1469    xdata = np.linspace(5.6,5.8,800)
1470    dx = xdata[1]-xdata[0]
1471    hplot = plotter.add('derivatives test for '+name).gca()
1472    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1473    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1474    myDict[name] += delt
1475    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1476    hplot.plot(xdata,(y1-y0)/delt,'r+')
1477
1478if __name__ == '__main__':
1479    import GSASIItestplot as plot
1480    global plotter
1481    plotter = plot.PlotNotebook()
1482#    test0()
1483#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1484    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1485        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1486        test2(name,shft)
1487    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1488        test3(name,shft)
1489    print "OK"
1490    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.