source: trunk/GSASIIpwd.py @ 1248

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

add setscale for SASD data
don't square the 1/cos(2-theta) correction to integrated intensities in ImageIntegrate?
scale SASD error bars by Scale
replace ':' with ';' in BkPk? parameter names
fix bug of missing Residuals in LS Refine

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 59.0 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII powder calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2014-03-13 19:16:20 +0000 (Thu, 13 Mar 2014) $
10# $Author: vondreele $
11# $Revision: 1248 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1248 2014-03-13 19:16:20Z 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: 1248 $")
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 getEpsVoigt(pos,alp,bet,sig,gam,xdata):
716    'needs a doc string'
717    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
718    Df /= np.sum(Df)
719    return Df 
720   
721def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
722    'needs a doc string'
723    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
724    sumDf = np.sum(Df)
725    return Df,dFdp,dFda,dFdb,dFds,dFdg   
726
727def ellipseSize(H,Sij,GB):
728    'needs a doc string'
729    HX = np.inner(H.T,GB)
730    lenHX = np.sqrt(np.sum(HX**2))
731    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
732    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
733    lenR = np.sqrt(np.sum(R**2))
734    return lenR
735
736def ellipseSizeDerv(H,Sij,GB):
737    'needs a doc string'
738    lenR = ellipseSize(H,Sij,GB)
739    delt = 0.001
740    dRdS = np.zeros(6)
741    for i in range(6):
742        dSij = Sij[:]
743        dSij[i] += delt
744        dRdS[i] = (ellipseSize(H,dSij,GB)-lenR)/delt
745    return lenR,dRdS
746
747def getHKLpeak(dmin,SGData,A):
748    'needs a doc string'
749    HKL = G2lat.GenHLaue(dmin,SGData,A)       
750    HKLs = []
751    for h,k,l,d in HKL:
752        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
753        if not ext:
754            HKLs.append([h,k,l,d,-1])
755    return HKLs
756
757def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
758    'needs a doc string'
759   
760    yb = getBackground('',parmDict,bakType,xdata)
761    yc = np.zeros_like(yb)
762    cw = np.diff(xdata)
763    cw = np.append(cw,cw[-1])
764    if 'C' in dataType:
765        shl = max(parmDict['SH/L'],0.002)
766        Ka2 = False
767        if 'Lam1' in parmDict.keys():
768            Ka2 = True
769            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
770            kRatio = parmDict['I(L2)/I(L1)']
771        iPeak = 0
772        while True:
773            try:
774                pos = parmDict['pos'+str(iPeak)]
775                theta = (pos-parmDict['Zero'])/2.0
776                intens = parmDict['int'+str(iPeak)]
777                sigName = 'sig'+str(iPeak)
778                if sigName in varyList:
779                    sig = parmDict[sigName]
780                else:
781                    sig = G2mth.getCWsig(parmDict,theta)
782                sig = max(sig,0.001)          #avoid neg sigma
783                gamName = 'gam'+str(iPeak)
784                if gamName in varyList:
785                    gam = parmDict[gamName]
786                else:
787                    gam = G2mth.getCWgam(parmDict,theta)
788                gam = max(gam,0.001)             #avoid neg gamma
789                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
790                iBeg = np.searchsorted(xdata,pos-fmin)
791                iFin = np.searchsorted(xdata,pos+fmin)
792                if not iBeg+iFin:       #peak below low limit
793                    iPeak += 1
794                    continue
795                elif not iBeg-iFin:     #peak above high limit
796                    return yb+yc
797                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
798                if Ka2:
799                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
800                    iBeg = np.searchsorted(xdata,pos2-fmin)
801                    iFin = np.searchsorted(xdata,pos2+fmin)
802                    if iBeg-iFin:
803                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
804                iPeak += 1
805            except KeyError:        #no more peaks to process
806                return yb+yc
807    else:
808        Pdabc = parmDict['Pdabc']
809        difC = parmDict['difC']
810        iPeak = 0
811        while True:
812            try:
813                pos = parmDict['pos'+str(iPeak)]               
814                tof = pos-parmDict['Zero']
815                dsp = tof/difC
816                intens = parmDict['int'+str(iPeak)]
817                alpName = 'alp'+str(iPeak)
818                if alpName in varyList:
819                    alp = parmDict[alpName]
820                else:
821                    if len(Pdabc):
822                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
823                    else:
824                        alp = G2mth.getTOFalpha(parmDict,dsp)
825                betName = 'bet'+str(iPeak)
826                if betName in varyList:
827                    bet = parmDict[betName]
828                else:
829                    if len(Pdabc):
830                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
831                    else:
832                        bet = G2mth.getTOFbeta(parmDict,dsp)
833                sigName = 'sig'+str(iPeak)
834                if sigName in varyList:
835                    sig = parmDict[sigName]
836                else:
837                    sig = G2mth.getTOFsig(parmDict,dsp)
838                gamName = 'gam'+str(iPeak)
839                if gamName in varyList:
840                    gam = parmDict[gamName]
841                else:
842                    gam = G2mth.getTOFgamma(parmDict,dsp)
843                gam = max(gam,0.001)             #avoid neg gamma
844                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
845                iBeg = np.searchsorted(xdata,pos-fmin)
846                iFin = np.searchsorted(xdata,pos+fmax)
847                lenX = len(xdata)
848                if not iBeg:
849                    iFin = np.searchsorted(xdata,pos+fmax)
850                elif iBeg == lenX:
851                    iFin = iBeg
852                else:
853                    iFin = np.searchsorted(xdata,pos+fmax)
854                if not iBeg+iFin:       #peak below low limit
855                    iPeak += 1
856                    continue
857                elif not iBeg-iFin:     #peak above high limit
858                    return yb+yc
859                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
860                iPeak += 1
861            except KeyError:        #no more peaks to process
862                return yb+yc
863           
864def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
865    'needs a doc string'
866# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
867    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
868    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
869    if 'Back:0' in varyList:            #background derivs are in front if present
870        dMdv[0:len(dMdb)] = dMdb
871    names = ['DebyeA','DebyeR','DebyeU']
872    for name in varyList:
873        if 'Debye' in name:
874            parm,id = name.split(':')
875            ip = names.index(parm)
876            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
877    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
878    for name in varyList:
879        if 'BkPk' in name:
880            parm,id = name.split(':')
881            ip = names.index(parm)
882            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
883    cw = np.diff(xdata)
884    cw = np.append(cw,cw[-1])
885    if 'C' in dataType:
886        shl = max(parmDict['SH/L'],0.002)
887        Ka2 = False
888        if 'Lam1' in parmDict.keys():
889            Ka2 = True
890            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
891            kRatio = parmDict['I(L2)/I(L1)']
892        iPeak = 0
893        while True:
894            try:
895                pos = parmDict['pos'+str(iPeak)]
896                theta = (pos-parmDict['Zero'])/2.0
897                intens = parmDict['int'+str(iPeak)]
898                sigName = 'sig'+str(iPeak)
899                tanth = tand(theta)
900                costh = cosd(theta)
901                if sigName in varyList:
902                    sig = parmDict[sigName]
903                else:
904                    sig = G2mth.getCWsig(parmDict,theta)
905                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(theta)
906                sig = max(sig,0.001)          #avoid neg sigma
907                gamName = 'gam'+str(iPeak)
908                if gamName in varyList:
909                    gam = parmDict[gamName]
910                else:
911                    gam = G2mth.getCWgam(parmDict,theta)
912                    dgdX,dgdY = G2mth.getCWgamDeriv(theta)
913                gam = max(gam,0.001)             #avoid neg gamma
914                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
915                iBeg = np.searchsorted(xdata,pos-fmin)
916                iFin = np.searchsorted(xdata,pos+fmin)
917                if not iBeg+iFin:       #peak below low limit
918                    iPeak += 1
919                    continue
920                elif not iBeg-iFin:     #peak above high limit
921                    break
922                dMdpk = np.zeros(shape=(6,len(xdata)))
923                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
924                for i in range(1,5):
925                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
926                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
927                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
928                if Ka2:
929                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
930                    iBeg = np.searchsorted(xdata,pos2-fmin)
931                    iFin = np.searchsorted(xdata,pos2+fmin)
932                    if iBeg-iFin:
933                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
934                        for i in range(1,5):
935                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
936                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
937                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
938                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
939                for parmName in ['pos','int','sig','gam']:
940                    try:
941                        idx = varyList.index(parmName+str(iPeak))
942                        dMdv[idx] = dervDict[parmName]
943                    except ValueError:
944                        pass
945                if 'U' in varyList:
946                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
947                if 'V' in varyList:
948                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
949                if 'W' in varyList:
950                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
951                if 'X' in varyList:
952                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
953                if 'Y' in varyList:
954                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
955                if 'SH/L' in varyList:
956                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
957                if 'I(L2)/I(L1)' in varyList:
958                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
959                iPeak += 1
960            except KeyError:        #no more peaks to process
961                break
962    else:
963        Pdabc = parmDict['Pdabc']
964        difC = parmDict['difC']
965        iPeak = 0
966        while True:
967            try:
968                pos = parmDict['pos'+str(iPeak)]               
969                tof = pos-parmDict['Zero']
970                dsp = tof/difC
971                intens = parmDict['int'+str(iPeak)]
972                alpName = 'alp'+str(iPeak)
973                if alpName in varyList:
974                    alp = parmDict[alpName]
975                else:
976                    if len(Pdabc):
977                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
978                    else:
979                        alp = G2mth.getTOFalpha(parmDict,dsp)
980                        dada0 = G2mth.getTOFalphaDeriv(dsp)
981                betName = 'bet'+str(iPeak)
982                if betName in varyList:
983                    bet = parmDict[betName]
984                else:
985                    if len(Pdabc):
986                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
987                    else:
988                        bet = G2mth.getTOFbeta(parmDict,dsp)
989                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
990                sigName = 'sig'+str(iPeak)
991                if sigName in varyList:
992                    sig = parmDict[sigName]
993                else:
994                    sig = G2mth.getTOFsig(parmDict,dsp)
995                    dsds0,dsds1,dsds2 = G2mth.getTOFsigDeriv(dsp)
996                gamName = 'gam'+str(iPeak)
997                if gamName in varyList:
998                    gam = parmDict[gamName]
999                else:
1000                    gam = G2mth.getTOFgamma(parmDict,dsp)
1001                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1002                gam = max(gam,0.001)             #avoid neg gamma
1003                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1004                iBeg = np.searchsorted(xdata,pos-fmin)
1005                lenX = len(xdata)
1006                if not iBeg:
1007                    iFin = np.searchsorted(xdata,pos+fmax)
1008                elif iBeg == lenX:
1009                    iFin = iBeg
1010                else:
1011                    iFin = np.searchsorted(xdata,pos+fmax)
1012                if not iBeg+iFin:       #peak below low limit
1013                    iPeak += 1
1014                    continue
1015                elif not iBeg-iFin:     #peak above high limit
1016                    break
1017                dMdpk = np.zeros(shape=(7,len(xdata)))
1018                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1019                for i in range(1,6):
1020                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1021                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1022                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1023                for parmName in ['pos','int','alp','bet','sig','gam']:
1024                    try:
1025                        idx = varyList.index(parmName+str(iPeak))
1026                        dMdv[idx] = dervDict[parmName]
1027                    except ValueError:
1028                        pass
1029                if 'alpha' in varyList:
1030                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1031                if 'beta-0' in varyList:
1032                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1033                if 'beta-1' in varyList:
1034                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1035                if 'beta-q' in varyList:
1036                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1037                if 'sig-0' in varyList:
1038                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1039                if 'sig-1' in varyList:
1040                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1041                if 'sig-q' in varyList:
1042                    dMdv[varyList.index('sig-q')] += dsds2*dervDict['sig']
1043                if 'X' in varyList:
1044                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1045                if 'Y' in varyList:
1046                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1047                iPeak += 1
1048            except KeyError:        #no more peaks to process
1049                break
1050    return dMdv
1051       
1052def Dict2Values(parmdict, varylist):
1053    '''Use before call to leastsq to setup list of values for the parameters
1054    in parmdict, as selected by key in varylist'''
1055    return [parmdict[key] for key in varylist] 
1056   
1057def Values2Dict(parmdict, varylist, values):
1058    ''' Use after call to leastsq to update the parameter dictionary with
1059    values corresponding to keys in varylist'''
1060    parmdict.update(zip(varylist,values))
1061   
1062def SetBackgroundParms(Background):
1063    'needs a doc string'
1064    if len(Background) == 1:            # fix up old backgrounds
1065        BackGround.append({'nDebye':0,'debyeTerms':[]})
1066    bakType,bakFlag = Background[0][:2]
1067    backVals = Background[0][3:]
1068    backNames = ['Back:'+str(i) for i in range(len(backVals))]
1069    Debye = Background[1]           #also has background peaks stuff
1070    backDict = dict(zip(backNames,backVals))
1071    backVary = []
1072    if bakFlag:
1073        backVary = backNames
1074
1075    backDict['nDebye'] = Debye['nDebye']
1076    debyeDict = {}
1077    debyeList = []
1078    for i in range(Debye['nDebye']):
1079        debyeNames = ['DebyeA:'+str(i),'DebyeR:'+str(i),'DebyeU:'+str(i)]
1080        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1081        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1082    debyeVary = []
1083    for item in debyeList:
1084        if item[1]:
1085            debyeVary.append(item[0])
1086    backDict.update(debyeDict)
1087    backVary += debyeVary
1088
1089    backDict['nPeaks'] = Debye['nPeaks']
1090    peaksDict = {}
1091    peaksList = []
1092    for i in range(Debye['nPeaks']):
1093        peaksNames = ['BkPkpos:'+str(i),'BkPkint:'+str(i),'BkPksig:'+str(i),'BkPkgam:'+str(i)]
1094        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1095        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1096    peaksVary = []
1097    for item in peaksList:
1098        if item[1]:
1099            peaksVary.append(item[0])
1100    backDict.update(peaksDict)
1101    backVary += peaksVary   
1102    return bakType,backDict,backVary
1103           
1104def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
1105    'needs a doc string'
1106   
1107       
1108    def GetBackgroundParms(parmList,Background):
1109        iBak = 0
1110        while True:
1111            try:
1112                bakName = 'Back:'+str(iBak)
1113                Background[0][iBak+3] = parmList[bakName]
1114                iBak += 1
1115            except KeyError:
1116                break
1117        iDb = 0
1118        while True:
1119            names = ['DebyeA:','DebyeR:','DebyeU:']
1120            try:
1121                for i,name in enumerate(names):
1122                    val = parmList[name+str(iDb)]
1123                    Background[1]['debyeTerms'][iDb][2*i] = val
1124                iDb += 1
1125            except KeyError:
1126                break
1127        iDb = 0
1128        while True:
1129            names = ['BkPkpos:','BkPkint:','BkPksig:','BkPkgam:']
1130            try:
1131                for i,name in enumerate(names):
1132                    val = parmList[name+str(iDb)]
1133                    Background[1]['peaksList'][iDb][2*i] = val
1134                iDb += 1
1135            except KeyError:
1136                break
1137               
1138    def BackgroundPrint(Background,sigDict):
1139        if Background[0][1]:
1140            print 'Background coefficients for',Background[0][0],'function'
1141            ptfmt = "%12.5f"
1142            ptstr =  'values:'
1143            sigstr = 'esds  :'
1144            for i,back in enumerate(Background[0][3:]):
1145                ptstr += ptfmt % (back)
1146                sigstr += ptfmt % (sigDict['Back:'+str(i)])
1147            print ptstr
1148            print sigstr
1149        else:
1150            print 'Background not refined'
1151        if Background[1]['nDebye']:
1152            parms = ['DebyeA','DebyeR','DebyeU']
1153            print 'Debye diffuse scattering coefficients'
1154            ptfmt = "%12.5f"
1155            names =   'names :'
1156            ptstr =  'values:'
1157            sigstr = 'esds  :'
1158            for item in sigDict:
1159                if 'Debye' in item:
1160                    names += '%12s'%(item)
1161                    sigstr += ptfmt%(sigDict[item])
1162                    parm,id = item.split(':')
1163                    ip = parms.index(parm)
1164                    ptstr += ptfmt%(Background[1]['debyeTerms'][int(id)][2*ip])
1165            print names
1166            print ptstr
1167            print sigstr
1168        if Background[1]['nPeaks']:
1169            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1170            print 'Peaks in background coefficients'
1171            ptfmt = "%12.5f"
1172            names =   'names :'
1173            ptstr =  'values:'
1174            sigstr = 'esds  :'
1175            for item in sigDict:
1176                if 'BkPk' in item:
1177                    names += '%12s'%(item)
1178                    sigstr += ptfmt%(sigDict[item])
1179                    parm,id = item.split(':')
1180                    ip = parms.index(parm)
1181                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1182            print names
1183            print ptstr
1184            print sigstr
1185                           
1186    def SetInstParms(Inst):
1187        dataType = Inst['Type'][0]
1188        insVary = []
1189        insNames = []
1190        insVals = []
1191        for parm in Inst:
1192            insNames.append(parm)
1193            insVals.append(Inst[parm][1])
1194            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1195                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',] and Inst[parm][2]:
1196                    insVary.append(parm)
1197        instDict = dict(zip(insNames,insVals))
1198        instDict['X'] = max(instDict['X'],0.01)
1199        instDict['Y'] = max(instDict['Y'],0.01)
1200        if 'SH/L' in instDict:
1201            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1202        return dataType,instDict,insVary
1203       
1204    def GetInstParms(parmDict,Inst,varyList,Peaks):
1205        for name in Inst:
1206            Inst[name][1] = parmDict[name]
1207        iPeak = 0
1208        while True:
1209            try:
1210                sigName = 'sig'+str(iPeak)
1211                pos = parmDict['pos'+str(iPeak)]
1212                if sigName not in varyList:
1213                    if 'C' in Inst['Type'][0]:
1214                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1215                    else:
1216                        dsp = pos/Inst['difC'][1]
1217                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1218                gamName = 'gam'+str(iPeak)
1219                if gamName not in varyList:
1220                    if 'C' in Inst['Type'][0]:
1221                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1222                    else:
1223                        dsp = pos/Inst['difC'][1]
1224                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1225                iPeak += 1
1226            except KeyError:
1227                break
1228       
1229    def InstPrint(Inst,sigDict):
1230        print 'Instrument Parameters:'
1231        ptfmt = "%12.6f"
1232        ptlbls = 'names :'
1233        ptstr =  'values:'
1234        sigstr = 'esds  :'
1235        for parm in Inst:
1236            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1237                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-q',]:
1238                ptlbls += "%s" % (parm.center(12))
1239                ptstr += ptfmt % (Inst[parm][1])
1240                if parm in sigDict:
1241                    sigstr += ptfmt % (sigDict[parm])
1242                else:
1243                    sigstr += 12*' '
1244        print ptlbls
1245        print ptstr
1246        print sigstr
1247
1248    def SetPeaksParms(dataType,Peaks):
1249        peakNames = []
1250        peakVary = []
1251        peakVals = []
1252        if 'C' in dataType:
1253            names = ['pos','int','sig','gam']
1254        else:
1255            names = ['pos','int','alp','bet','sig','gam']
1256        for i,peak in enumerate(Peaks):
1257            for j,name in enumerate(names):
1258                peakVals.append(peak[2*j])
1259                parName = name+str(i)
1260                peakNames.append(parName)
1261                if peak[2*j+1]:
1262                    peakVary.append(parName)
1263        return dict(zip(peakNames,peakVals)),peakVary
1264               
1265    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1266        if 'C' in Inst['Type'][0]:
1267            names = ['pos','int','sig','gam']
1268        else:
1269            names = ['pos','int','alp','bet','sig','gam']
1270        for i,peak in enumerate(Peaks):
1271            pos = parmDict['pos'+str(i)]
1272            if 'difC' in Inst:
1273                dsp = pos/Inst['difC'][1]
1274            for j in range(len(names)):
1275                parName = names[j]+str(i)
1276                if parName in varyList:
1277                    peak[2*j] = parmDict[parName]
1278                elif 'alpha' in parName:
1279                    peak[2*j] = parmDict['alpha']/dsp
1280                elif 'beta' in parName:
1281                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1282                elif 'sig' in parName:
1283                    if 'C' in Inst['Type'][0]:
1284                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1285                    else:
1286                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1287                elif 'gam' in parName:
1288                    if 'C' in Inst['Type'][0]:
1289                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1290                    else:
1291                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1292                       
1293    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1294        print 'Peak coefficients:'
1295        if 'C' in dataType:
1296            names = ['pos','int','sig','gam']
1297        else:
1298            names = ['pos','int','alp','bet','sig','gam']           
1299        head = 13*' '
1300        for name in names:
1301            head += name.center(10)+'esd'.center(10)
1302        print head
1303        if 'C' in dataType:
1304            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1305        else:
1306            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%10.3f",'bet':"%10.5f",'sig':"%10.3f",'gam':"%10.3f"}
1307        for i,peak in enumerate(Peaks):
1308            ptstr =  ':'
1309            for j in range(len(names)):
1310                name = names[j]
1311                parName = name+str(i)
1312                ptstr += ptfmt[name] % (parmDict[parName])
1313                if parName in varyList:
1314#                    ptstr += G2IO.ValEsd(parmDict[parName],sigDict[parName])
1315                    ptstr += ptfmt[name] % (sigDict[parName])
1316                else:
1317#                    ptstr += G2IO.ValEsd(parmDict[parName],0.0)
1318                    ptstr += 10*' '
1319            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1320               
1321    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1322        parmdict.update(zip(varylist,values))
1323        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1324           
1325    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1326        parmdict.update(zip(varylist,values))
1327        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1328        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1329        if dlg:
1330            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1331            if not GoOn:
1332                return -M           #abort!!
1333        return M
1334       
1335    if controls:
1336        Ftol = controls['min dM/M']
1337        derivType = controls['deriv type']
1338    else:
1339        Ftol = 0.0001
1340        derivType = 'analytic'
1341    if oneCycle:
1342        Ftol = 1.0
1343    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1344    yc *= 0.                            #set calcd ones to zero
1345    yb *= 0.
1346    yd *= 0.
1347    xBeg = np.searchsorted(x,Limits[0])
1348    xFin = np.searchsorted(x,Limits[1])
1349    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1350    dataType,insDict,insVary = SetInstParms(Inst)
1351    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1352    parmDict = {}
1353    parmDict.update(bakDict)
1354    parmDict.update(insDict)
1355    parmDict.update(peakDict)
1356    parmDict['Pdabc'] = []      #dummy Pdabc
1357    parmDict.update(Inst2)      #put in real one if there
1358    varyList = bakVary+insVary+peakVary
1359    while True:
1360        begin = time.time()
1361        values =  np.array(Dict2Values(parmDict, varyList))
1362        if FitPgm == 'LSQ':
1363            try:
1364                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1365                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1366                ncyc = int(result[2]['nfev']/2)
1367            finally:
1368                dlg.Destroy()
1369            runtime = time.time()-begin   
1370            chisq = np.sum(result[2]['fvec']**2)
1371            Values2Dict(parmDict, varyList, result[0])
1372            Rwp = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1373            GOF = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1374            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1375            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1376            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1377            try:
1378                sig = np.sqrt(np.diag(result[1])*GOF)
1379                if np.any(np.isnan(sig)):
1380                    print '*** Least squares aborted - some invalid esds possible ***'
1381                break                   #refinement succeeded - finish up!
1382            except ValueError:          #result[1] is None on singular matrix
1383                print '**** Refinement failed - singular matrix ****'
1384                Ipvt = result[2]['ipvt']
1385                for i,ipvt in enumerate(Ipvt):
1386                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1387                        print 'Removing parameter: ',varyList[ipvt-1]
1388                        del(varyList[ipvt-1])
1389                        break
1390        elif FitPgm == 'BFGS':
1391            print 'Other program here'
1392            return
1393       
1394    sigDict = dict(zip(varyList,sig))
1395    yb[xBeg:xFin] = getBackground('',parmDict,bakType,x[xBeg:xFin])
1396    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1397    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1398    GetBackgroundParms(parmDict,Background)
1399    BackgroundPrint(Background,sigDict)
1400    GetInstParms(parmDict,Inst,varyList,Peaks)
1401    InstPrint(Inst,sigDict)
1402    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1403    PeaksPrint(dataType,parmDict,sigDict,varyList)
1404
1405def calcIncident(Iparm,xdata):
1406    'needs a doc string'
1407
1408    def IfunAdv(Iparm,xdata):
1409        Itype = Iparm['Itype']
1410        Icoef = Iparm['Icoeff']
1411        DYI = np.ones((12,xdata.shape[0]))
1412        YI = np.ones_like(xdata)*Icoef[0]
1413       
1414        x = xdata/1000.                 #expressions are in ms
1415        if Itype == 'Exponential':
1416            for i in range(1,10,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        elif 'Maxwell'in Itype:
1422            Eterm = np.exp(-Icoef[2]/x**2)
1423            DYI[1] = Eterm/x**5
1424            DYI[2] = -Icoef[1]*DYI[1]/x**2
1425            YI += (Icoef[1]*Eterm/x**5)
1426            if 'Exponential' in Itype:
1427                for i in range(3,12,2):
1428                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1429                    YI += Icoef[i]*Eterm
1430                    DYI[i] *= Eterm
1431                    DYI[i+1] *= -Icoef[i]*x**((i+1)/2)
1432            else:   #Chebyschev
1433                T = (2./x)-1.
1434                Ccof = np.ones((12,xdata.shape[0]))
1435                Ccof[1] = T
1436                for i in range(2,12):
1437                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1438                for i in range(1,10):
1439                    YI += Ccof[i]*Icoef[i+2]
1440                    DYI[i+2] =Ccof[i]
1441        return YI,DYI
1442       
1443    Iesd = np.array(Iparm['Iesd'])
1444    Icovar = Iparm['Icovar']
1445    YI,DYI = IfunAdv(Iparm,xdata)
1446    YI = np.where(YI>0,YI,1.)
1447    WYI = np.zeros_like(xdata)
1448    vcov = np.zeros((12,12))
1449    k = 0
1450    for i in range(12):
1451        for j in range(i,12):
1452            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1453            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1454            k += 1
1455    M = np.inner(vcov,DYI.T)
1456    WYI = np.sum(M*DYI,axis=0)
1457    WYI = np.where(WYI>0.,WYI,0.)
1458    return YI,WYI
1459   
1460#testing data
1461NeedTestData = True
1462def TestData():
1463    'needs a doc string'
1464#    global NeedTestData
1465    NeedTestData = False
1466    global bakType
1467    bakType = 'chebyschev'
1468    global xdata
1469    xdata = np.linspace(4.0,40.0,36000)
1470    global parmDict0
1471    parmDict0 = {
1472        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1473        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1474        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1475        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1476        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1477        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1478        }
1479    global parmDict1
1480    parmDict1 = {
1481        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1482        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1483        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1484        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1485        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1486        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1487        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1488        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1489        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1490        }
1491    global parmDict2
1492    parmDict2 = {
1493        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1494        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1495        'Back0':5.,'Back1':-0.02,'Back2':.004,
1496#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1497        }
1498    global varyList
1499    varyList = []
1500
1501def test0():
1502    if NeedTestData: TestData()
1503    msg = 'test '
1504    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1505    gplot.plot(xdata,getBackground('',parmDict0,bakType,xdata))   
1506    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1507    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1508    fplot.plot(xdata,getBackground('',parmDict1,bakType,xdata))   
1509    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1510   
1511def test1():
1512    if NeedTestData: TestData()
1513    time0 = time.time()
1514    for i in range(100):
1515        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1516    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1517   
1518def test2(name,delt):
1519    if NeedTestData: TestData()
1520    varyList = [name,]
1521    xdata = np.linspace(5.6,5.8,400)
1522    hplot = plotter.add('derivatives test for '+name).gca()
1523    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1524    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1525    parmDict2[name] += delt
1526    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1527    hplot.plot(xdata,(y1-y0)/delt,'r+')
1528   
1529def test3(name,delt):
1530    if NeedTestData: TestData()
1531    names = ['pos','sig','gam','shl']
1532    idx = names.index(name)
1533    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1534    xdata = np.linspace(5.6,5.8,800)
1535    dx = xdata[1]-xdata[0]
1536    hplot = plotter.add('derivatives test for '+name).gca()
1537    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1538    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1539    myDict[name] += delt
1540    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1541    hplot.plot(xdata,(y1-y0)/delt,'r+')
1542
1543if __name__ == '__main__':
1544    import GSASIItestplot as plot
1545    global plotter
1546    plotter = plot.PlotNotebook()
1547#    test0()
1548#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1549    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1550        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1551        test2(name,shft)
1552    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1553        test3(name,shft)
1554    print "OK"
1555    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.