source: trunk/GSASIIpwd.py @ 1176

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

add a bit of Doc to surface roughness
remove deselect tree line - fails in Windows

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