source: trunk/GSASIIpwd.py @ 1175

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

add surface roughness (Surotti model) to Bragg-Brentano sample parameters
seems to work OK but I don't have a good test data set.

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