source: trunk/GSASIIpwd.py @ 1373

Last change on this file since 1373 was 1373, checked in by toby, 8 years ago

minor source documention fixes

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