source: trunk/GSASIIpwd.py @ 1662

Last change on this file since 1662 was 1662, checked in by vondreele, 8 years ago

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