source: trunk/GSASIIpwd.py @ 808

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

Add Babinet modification to form factors - works OK
Use atom Ids for restraints
add 'all' in atom constraint editing routines
New add routines for bond restraints (angle add not done yet)
FindAtomIndexByIDs, FillAtomLookUp?, GetAtomsByID, & GetAtomItemsById? all now in GSASIImath.py
removed from GSASIIphsGUI.py
change phase tick mark colors
get GSASIIstruct.py to recognize protein atoms & do a protein calc.

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