source: trunk/GSASIIpwd.py @ 1474

Last change on this file since 1474 was 1474, checked in by vondreele, 7 years ago

1st MC/SA tutorial
various MC/SA fixes
fix to background peak fitting for CW & TOF

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