source: trunk/GSASIIpwd.py @ 1274

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

implement Bragg peaks in SASD - required new fortran routines
implement unified (Guinier + Porod) and Porod models in SASD
implement monodisperse models in SASD
correct bin width issue in lognormal, etc. fitting

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