source: trunk/GSASIIpwd.py @ 1412

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

fix theta vs 2-theta error in peak fitting; gave wrong sig,gam vs U,V,W,X,Y now all match OK

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 59.7 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-07-06 14:21:55 +0000 (Sun, 06 Jul 2014) $
10# $Author: vondreele $
11# $Revision: 1412 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1412 2014-07-06 14:21:55Z 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: 1412 $")
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                tth = (pos-parmDict['Zero'])
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,tth)
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,tth)
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                tth = (pos-parmDict['Zero'])
915                intens = parmDict['int'+str(iPeak)]
916                sigName = 'sig'+str(iPeak)
917                if sigName in varyList:
918                    sig = parmDict[sigName]
919                    dsdU = dsdV = dsdW = 0
920                else:
921                    sig = G2mth.getCWsig(parmDict,tth)
922                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
923                sig = max(sig,0.001)          #avoid neg sigma
924                gamName = 'gam'+str(iPeak)
925                if gamName in varyList:
926                    gam = parmDict[gamName]
927                    dgdX = dgdY = 0
928                else:
929                    gam = G2mth.getCWgam(parmDict,tth)
930                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
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                        dad0 = 0
997                    else:
998                        alp = G2mth.getTOFalpha(parmDict,dsp)
999                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1000                betName = 'bet'+str(iPeak)
1001                if betName in varyList:
1002                    bet = parmDict[betName]
1003                else:
1004                    if len(Pdabc):
1005                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1006                        dbdb0 = dbdb1 = dbdb2 = 0
1007                    else:
1008                        bet = G2mth.getTOFbeta(parmDict,dsp)
1009                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1010                sigName = 'sig'+str(iPeak)
1011                if sigName in varyList:
1012                    sig = parmDict[sigName]
1013                    dsds0 = dsds1 = dsds2 = 0
1014                else:
1015                    sig = G2mth.getTOFsig(parmDict,dsp)
1016                    dsds0,dsds1,dsds2 = G2mth.getTOFsigDeriv(dsp)
1017                gamName = 'gam'+str(iPeak)
1018                if gamName in varyList:
1019                    gam = parmDict[gamName]
1020                    dsdX = dsdY = 0
1021                else:
1022                    gam = G2mth.getTOFgamma(parmDict,dsp)
1023                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1024                gam = max(gam,0.001)             #avoid neg gamma
1025                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1026                iBeg = np.searchsorted(xdata,pos-fmin)
1027                lenX = len(xdata)
1028                if not iBeg:
1029                    iFin = np.searchsorted(xdata,pos+fmax)
1030                elif iBeg == lenX:
1031                    iFin = iBeg
1032                else:
1033                    iFin = np.searchsorted(xdata,pos+fmax)
1034                if not iBeg+iFin:       #peak below low limit
1035                    iPeak += 1
1036                    continue
1037                elif not iBeg-iFin:     #peak above high limit
1038                    break
1039                dMdpk = np.zeros(shape=(7,len(xdata)))
1040                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1041                for i in range(1,6):
1042                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1043                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1044                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1045                for parmName in ['pos','int','alp','bet','sig','gam']:
1046                    try:
1047                        idx = varyList.index(parmName+str(iPeak))
1048                        dMdv[idx] = dervDict[parmName]
1049                    except ValueError:
1050                        pass
1051                if 'alpha' in varyList:
1052                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1053                if 'beta-0' in varyList:
1054                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1055                if 'beta-1' in varyList:
1056                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1057                if 'beta-q' in varyList:
1058                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1059                if 'sig-0' in varyList:
1060                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1061                if 'sig-1' in varyList:
1062                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1063                if 'sig-q' in varyList:
1064                    dMdv[varyList.index('sig-q')] += dsds2*dervDict['sig']
1065                if 'X' in varyList:
1066                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1067                if 'Y' in varyList:
1068                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1069                iPeak += 1
1070            except KeyError:        #no more peaks to process
1071                break
1072    return dMdv
1073       
1074def Dict2Values(parmdict, varylist):
1075    '''Use before call to leastsq to setup list of values for the parameters
1076    in parmdict, as selected by key in varylist'''
1077    return [parmdict[key] for key in varylist] 
1078   
1079def Values2Dict(parmdict, varylist, values):
1080    ''' Use after call to leastsq to update the parameter dictionary with
1081    values corresponding to keys in varylist'''
1082    parmdict.update(zip(varylist,values))
1083   
1084def SetBackgroundParms(Background):
1085    'needs a doc string'
1086    if len(Background) == 1:            # fix up old backgrounds
1087        BackGround.append({'nDebye':0,'debyeTerms':[]})
1088    bakType,bakFlag = Background[0][:2]
1089    backVals = Background[0][3:]
1090    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1091    Debye = Background[1]           #also has background peaks stuff
1092    backDict = dict(zip(backNames,backVals))
1093    backVary = []
1094    if bakFlag:
1095        backVary = backNames
1096
1097    backDict['nDebye'] = Debye['nDebye']
1098    debyeDict = {}
1099    debyeList = []
1100    for i in range(Debye['nDebye']):
1101        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1102        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1103        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1104    debyeVary = []
1105    for item in debyeList:
1106        if item[1]:
1107            debyeVary.append(item[0])
1108    backDict.update(debyeDict)
1109    backVary += debyeVary
1110
1111    backDict['nPeaks'] = Debye['nPeaks']
1112    peaksDict = {}
1113    peaksList = []
1114    for i in range(Debye['nPeaks']):
1115        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1116        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1117        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1118    peaksVary = []
1119    for item in peaksList:
1120        if item[1]:
1121            peaksVary.append(item[0])
1122    backDict.update(peaksDict)
1123    backVary += peaksVary   
1124    return bakType,backDict,backVary
1125           
1126def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1127    'needs a doc string'
1128       
1129    def GetBackgroundParms(parmList,Background):
1130        iBak = 0
1131        while True:
1132            try:
1133                bakName = 'Back:'+str(iBak)
1134                Background[0][iBak+3] = parmList[bakName]
1135                iBak += 1
1136            except KeyError:
1137                break
1138        iDb = 0
1139        while True:
1140            names = ['DebyeA:','DebyeR:','DebyeU:']
1141            try:
1142                for i,name in enumerate(names):
1143                    val = parmList[name+str(iDb)]
1144                    Background[1]['debyeTerms'][iDb][2*i] = val
1145                iDb += 1
1146            except KeyError:
1147                break
1148        iDb = 0
1149        while True:
1150            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1151            try:
1152                for i,name in enumerate(names):
1153                    val = parmList[name+str(iDb)]
1154                    Background[1]['peaksList'][iDb][2*i] = val
1155                iDb += 1
1156            except KeyError:
1157                break
1158               
1159    def BackgroundPrint(Background,sigDict):
1160        if Background[0][1]:
1161            print 'Background coefficients for',Background[0][0],'function'
1162            ptfmt = "%12.5f"
1163            ptstr =  'values:'
1164            sigstr = 'esds  :'
1165            for i,back in enumerate(Background[0][3:]):
1166                ptstr += ptfmt % (back)
1167                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1168            print ptstr
1169            print sigstr
1170        else:
1171            print 'Background not refined'
1172        if Background[1]['nDebye']:
1173            parms = ['DebyeA','DebyeR','DebyeU']
1174            print 'Debye diffuse scattering coefficients'
1175            ptfmt = "%12.5f"
1176            names =   'names :'
1177            ptstr =  'values:'
1178            sigstr = 'esds  :'
1179            for item in sigDict:
1180                if 'Debye' in item:
1181                    names += '%12s'%(item)
1182                    sigstr += ptfmt%(sigDict[item])
1183                    parm,id = item.split(':')
1184                    ip = parms.index(parm)
1185                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1186            print names
1187            print ptstr
1188            print sigstr
1189        if Background[1]['nPeaks']:
1190            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1191            print 'Peaks in background coefficients'
1192            ptfmt = "%12.5f"
1193            names =   'names :'
1194            ptstr =  'values:'
1195            sigstr = 'esds  :'
1196            for item in sigDict:
1197                if 'BkPk' in item:
1198                    names += '%12s'%(item)
1199                    sigstr += ptfmt%(sigDict[item])
1200                    parm,id = item.split(':')
1201                    ip = parms.index(parm)
1202                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1203            print names
1204            print ptstr
1205            print sigstr
1206                           
1207    def SetInstParms(Inst):
1208        dataType = Inst['Type'][0]
1209        insVary = []
1210        insNames = []
1211        insVals = []
1212        for parm in Inst:
1213            insNames.append(parm)
1214            insVals.append(Inst[parm][1])
1215            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1216                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',] and Inst[parm][2]:
1217                    insVary.append(parm)
1218        instDict = dict(zip(insNames,insVals))
1219        instDict['X'] = max(instDict['X'],0.01)
1220        instDict['Y'] = max(instDict['Y'],0.01)
1221        if 'SH/L' in instDict:
1222            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1223        return dataType,instDict,insVary
1224       
1225    def GetInstParms(parmDict,Inst,varyList,Peaks):
1226        for name in Inst:
1227            Inst[name][1] = parmDict[name]
1228        iPeak = 0
1229        while True:
1230            try:
1231                sigName = 'sig'+str(iPeak)
1232                pos = parmDict['pos'+str(iPeak)]
1233                if sigName not in varyList:
1234                    if 'C' in Inst['Type'][0]:
1235                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1236                    else:
1237                        dsp = pos/Inst['difC'][1]
1238                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1239                gamName = 'gam'+str(iPeak)
1240                if gamName not in varyList:
1241                    if 'C' in Inst['Type'][0]:
1242                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1243                    else:
1244                        dsp = pos/Inst['difC'][1]
1245                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1246                iPeak += 1
1247            except KeyError:
1248                break
1249       
1250    def InstPrint(Inst,sigDict):
1251        print 'Instrument Parameters:'
1252        ptfmt = "%12.6f"
1253        ptlbls = 'names :'
1254        ptstr =  'values:'
1255        sigstr = 'esds  :'
1256        for parm in Inst:
1257            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1258                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',]:
1259                ptlbls += "%s" % (parm.center(12))
1260                ptstr += ptfmt % (Inst[parm][1])
1261                if parm in sigDict:
1262                    sigstr += ptfmt % (sigDict[parm])
1263                else:
1264                    sigstr += 12*' '
1265        print ptlbls
1266        print ptstr
1267        print sigstr
1268
1269    def SetPeaksParms(dataType,Peaks):
1270        peakNames = []
1271        peakVary = []
1272        peakVals = []
1273        if 'C' in dataType:
1274            names = ['pos','int','sig','gam']
1275        else:
1276            names = ['pos','int','alp','bet','sig','gam']
1277        for i,peak in enumerate(Peaks):
1278            for j,name in enumerate(names):
1279                peakVals.append(peak[2*j])
1280                parName = name+str(i)
1281                peakNames.append(parName)
1282                if peak[2*j+1]:
1283                    peakVary.append(parName)
1284        return dict(zip(peakNames,peakVals)),peakVary
1285               
1286    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1287        if 'C' in Inst['Type'][0]:
1288            names = ['pos','int','sig','gam']
1289        else:   #'T'
1290            names = ['pos','int','alp','bet','sig','gam']
1291        for i,peak in enumerate(Peaks):
1292            pos = parmDict['pos'+str(i)]
1293            if 'difC' in Inst:
1294                dsp = pos/Inst['difC'][1]
1295            for j in range(len(names)):
1296                parName = names[j]+str(i)
1297                if parName in varyList:
1298                    peak[2*j] = parmDict[parName]
1299                elif 'alpha' in parName:
1300                    peak[2*j] = parmDict['alpha']/dsp
1301                elif 'beta' in parName:
1302                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1303                elif 'sig' in parName:
1304                    if 'C' in Inst['Type'][0]:
1305                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1306                    else:
1307                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1308                elif 'gam' in parName:
1309                    if 'C' in Inst['Type'][0]:
1310                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1311                    else:
1312                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1313                       
1314    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1315        print 'Peak coefficients:'
1316        if 'C' in dataType:
1317            names = ['pos','int','sig','gam']
1318        else:   #'T'
1319            names = ['pos','int','alp','bet','sig','gam']           
1320        head = 13*' '
1321        for name in names:
1322            head += name.center(10)+'esd'.center(10)
1323        print head
1324        if 'C' in dataType:
1325            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1326        else:
1327            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1328        for i,peak in enumerate(Peaks):
1329            ptstr =  ':'
1330            for j in range(len(names)):
1331                name = names[j]
1332                parName = name+str(i)
1333                ptstr += ptfmt[name] % (parmDict[parName])
1334                if parName in varyList:
1335#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1336                    ptstr += ptfmt[name] % (sigDict[parName])
1337                else:
1338#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1339                    ptstr += 10*' '
1340            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1341               
1342    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1343        parmdict.update(zip(varylist,values))
1344        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1345           
1346    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1347        parmdict.update(zip(varylist,values))
1348        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1349        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1350        if dlg:
1351            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1352            if not GoOn:
1353                return -M           #abort!!
1354        return M
1355       
1356    if controls:
1357        Ftol = controls['min dM/M']
1358        derivType = controls['deriv type']
1359    else:
1360        Ftol = 0.0001
1361        derivType = 'analytic'
1362    if oneCycle:
1363        Ftol = 1.0
1364    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1365    yc *= 0.                            #set calcd ones to zero
1366    yb *= 0.
1367    yd *= 0.
1368    xBeg = np.searchsorted(x,Limits[0])
1369    xFin = np.searchsorted(x,Limits[1])
1370    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1371    dataType,insDict,insVary = SetInstParms(Inst)
1372    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1373    parmDict = {}
1374    parmDict.update(bakDict)
1375    parmDict.update(insDict)
1376    parmDict.update(peakDict)
1377    parmDict['Pdabc'] = []      #dummy Pdabc
1378    parmDict.update(Inst2)      #put in real one if there
1379    varyList = bakVary+insVary+peakVary
1380    while True:
1381        begin = time.time()
1382        values =  np.array(Dict2Values(parmDict, varyList))
1383        if FitPgm == 'LSQ':
1384            try:
1385                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1386                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1387                ncyc = int(result[2]['nfev']/2)
1388            finally:
1389                dlg.Destroy()
1390            runtime = time.time()-begin   
1391            chisq = np.sum(result[2]['fvec']**2)
1392            Values2Dict(parmDict, varyList, result[0])
1393            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1394            GOF = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1395            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1396            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1397            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1398            try:
1399                sig = np.sqrt(np.diag(result[1])*GOF)
1400                if np.any(np.isnan(sig)):
1401                    print '*** Least squares aborted - some invalid esds possible ***'
1402                break                   #refinement succeeded - finish up!
1403            except ValueError:          #result[1] is None on singular matrix
1404                print '**** Refinement failed - singular matrix ****'
1405                Ipvt = result[2]['ipvt']
1406                for i,ipvt in enumerate(Ipvt):
1407                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1408                        print 'Removing parameter: ',varyList[ipvt-1]
1409                        del(varyList[ipvt-1])
1410                        break
1411        elif FitPgm == 'BFGS':
1412            print 'Other program here'
1413            return
1414       
1415    sigDict = dict(zip(varyList,sig))
1416    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1417    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1418    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1419    GetBackgroundParms(parmDict,Background)
1420    BackgroundPrint(Background,sigDict)
1421    GetInstParms(parmDict,Inst,varyList,Peaks)
1422    InstPrint(Inst,sigDict)
1423    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1424    PeaksPrint(dataType,parmDict,sigDict,varyList)
1425
1426def calcIncident(Iparm,xdata):
1427    'needs a doc string'
1428
1429    def IfunAdv(Iparm,xdata):
1430        Itype = Iparm['Itype']
1431        Icoef = Iparm['Icoeff']
1432        DYI = np.ones((12,xdata.shape[0]))
1433        YI = np.ones_like(xdata)*Icoef[0]
1434       
1435        x = xdata/1000.                 #expressions are in ms
1436        if Itype == 'Exponential':
1437            for i in range(1,10,2):
1438                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1439                YI += Icoef[i]*Eterm
1440                DYI[i] *= Eterm
1441                DYI[i+1] *= -Icoef[i]*x**((i+1)/2)           
1442        elif 'Maxwell'in Itype:
1443            Eterm = np.exp(-Icoef[2]/x**2)
1444            DYI[1] = Eterm/x**5
1445            DYI[2] = -Icoef[1]*DYI[1]/x**2
1446            YI += (Icoef[1]*Eterm/x**5)
1447            if 'Exponential' in Itype:
1448                for i in range(3,12,2):
1449                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1450                    YI += Icoef[i]*Eterm
1451                    DYI[i] *= Eterm
1452                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1453            else:   #Chebyschev
1454                T = (2./x)-1.
1455                Ccof = np.ones((12,xdata.shape[0]))
1456                Ccof[1] = T
1457                for i in range(2,12):
1458                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1459                for i in range(1,10):
1460                    YI += Ccof[i]*Icoef[i+2]
1461                    DYI[i+2] =Ccof[i]
1462        return YI,DYI
1463       
1464    Iesd = np.array(Iparm['Iesd'])
1465    Icovar = Iparm['Icovar']
1466    YI,DYI = IfunAdv(Iparm,xdata)
1467    YI = np.where(YI>0,YI,1.)
1468    WYI = np.zeros_like(xdata)
1469    vcov = np.zeros((12,12))
1470    k = 0
1471    for i in range(12):
1472        for j in range(i,12):
1473            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1474            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1475            k += 1
1476    M = np.inner(vcov,DYI.T)
1477    WYI = np.sum(M*DYI,axis=0)
1478    WYI = np.where(WYI>0.,WYI,0.)
1479    return YI,WYI
1480   
1481#testing data
1482NeedTestData = True
1483def TestData():
1484    'needs a doc string'
1485#    global NeedTestData
1486    NeedTestData = False
1487    global bakType
1488    bakType = 'chebyschev'
1489    global xdata
1490    xdata = np.linspace(4.0,40.0,36000)
1491    global parmDict0
1492    parmDict0 = {
1493        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1494        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1495        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1496        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1497        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1498        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1499        }
1500    global parmDict1
1501    parmDict1 = {
1502        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1503        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1504        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1505        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1506        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1507        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1508        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1509        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1510        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1511        }
1512    global parmDict2
1513    parmDict2 = {
1514        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1515        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1516        'Back0':5.,'Back1':-0.02,'Back2':.004,
1517#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1518        }
1519    global varyList
1520    varyList = []
1521
1522def test0():
1523    if NeedTestData: TestData()
1524    msg = 'test '
1525    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1526    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1527    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1528    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1529    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1530    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1531   
1532def test1():
1533    if NeedTestData: TestData()
1534    time0 = time.time()
1535    for i in range(100):
1536        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1537    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1538   
1539def test2(name,delt):
1540    if NeedTestData: TestData()
1541    varyList = [name,]
1542    xdata = np.linspace(5.6,5.8,400)
1543    hplot = plotter.add('derivatives test for '+name).gca()
1544    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1545    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1546    parmDict2[name] += delt
1547    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1548    hplot.plot(xdata,(y1-y0)/delt,'r+')
1549   
1550def test3(name,delt):
1551    if NeedTestData: TestData()
1552    names = ['pos','sig','gam','shl']
1553    idx = names.index(name)
1554    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1555    xdata = np.linspace(5.6,5.8,800)
1556    dx = xdata[1]-xdata[0]
1557    hplot = plotter.add('derivatives test for '+name).gca()
1558    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1559    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1560    myDict[name] += delt
1561    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1562    hplot.plot(xdata,(y1-y0)/delt,'r+')
1563
1564if __name__ == '__main__':
1565    import GSASIItestplot as plot
1566    global plotter
1567    plotter = plot.PlotNotebook()
1568#    test0()
1569#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1570    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1571        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1572        test2(name,shft)
1573    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1574        test3(name,shft)
1575    print "OK"
1576    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.