source: trunk/GSASIIpwd.py @ 997

Last change on this file since 997 was 997, checked in by vondreele, 10 years ago

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