source: trunk/GSASIIpwd.py @ 945

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

some mods for POWGEN data for single peak fitting
unify profile parm calcs into G2mth

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