source: trunk/GSASIIpwd.py @ 1462

Last change on this file since 1462 was 1462, checked in by vondreele, 9 years ago

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