source: trunk/GSASIIpwd.py @ 1277

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

1st working SASD fit version. Implement UnDo? option for SASD fitting.

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