source: trunk/GSASIIpwd.py @ 1878

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

refactor DDataGUI - mostly move event routines to be inside respective sizer routines
Enable Flack parameter - function OK; derivatives need work
Allow inversion of noncentrosymmetric structures in SymOpDialog? - to test enantiomers
some work to make pdf stuff neutron TOF friendly - not complete
fix Debye background function - now works; refactor result printing for it

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 67.6 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: 2015-06-04 18:49:58 +0000 (Thu, 04 Jun 2015) $
10# $Author: vondreele $
11# $Revision: 1878 $
12# $URL: trunk/GSASIIpwd.py $
13# $Id: GSASIIpwd.py 1878 2015-06-04 18:49:58Z 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: 1878 $")
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    if 'X' in inst['Type'][0]:
273        xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
274    if data['DetType'] == 'Image plate':
275        xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
276    XY = xydata['IofQ'][1]   
277    #convert to Q
278    hc = 12.397639
279    wave = G2mth.getWave(inst)
280    keV = hc/wave
281    minQ = npT2q(Tth[0],wave)
282    maxQ = npT2q(Tth[-1],wave)   
283    Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
284    dq = Qpoints[1]-Qpoints[0]
285    XY[0] = npT2q(XY[0],wave)   
286#    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
287    T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=0.0)      #OK for scipy 0.8
288    Qdata = T(Qpoints)
289   
290    qLimits = data['QScaleLim']
291    minQ = np.searchsorted(Qpoints,qLimits[0])
292    maxQ = np.searchsorted(Qpoints,qLimits[1])
293    newdata = []
294    xydata['IofQ'][1][0] = Qpoints
295    xydata['IofQ'][1][1] = Qdata
296    for item in xydata['IofQ'][1]:
297        newdata.append(item[:maxQ])
298    xydata['IofQ'][1] = newdata
299   
300
301    xydata['SofQ'] = copy.deepcopy(xydata['IofQ'])
302    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
303    Q = xydata['SofQ'][1][0]
304    ruland = Ruland(data['Ruland'],wave,Q,CF)
305#    auxPlot.append([Q,ruland,'Ruland'])     
306    CF *= ruland
307#    auxPlot.append([Q,CF,'CF-Corr'])
308    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
309    xydata['SofQ'][1][1] *= scale
310    xydata['SofQ'][1][1] -= CF
311    xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF
312    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
313    xydata['SofQ'][1][1] *= scale
314   
315    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
316    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
317    if data['Lorch']:
318        xydata['FofQ'][1][1] *= LorchWeight(Q)
319   
320    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
321    nR = len(xydata['GofR'][1][1])
322    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],4*nR)[:nR])
323    xydata['GofR'][1][0] = 0.5*np.pi*np.linspace(0,nR,nR)/qLimits[1]
324   
325       
326    return auxPlot
327       
328#GSASII peak fitting routines: Finger, Cox & Jephcoat model       
329
330def factorize(num):
331    ''' Provide prime number factors for integer num
332    :returns: dictionary of prime factors (keys) & power for each (data)
333    '''
334    factors = {}
335    orig = num
336
337    # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1
338    i, sqi = 2, 4
339    while sqi <= num:
340        while not num%i:
341            num /= i
342            factors[i] = factors.get(i, 0) + 1
343
344        sqi += 2*i + 1
345        i += 1
346
347    if num != 1 and num != orig:
348        factors[num] = factors.get(num, 0) + 1
349
350    if factors:
351        return factors
352    else:
353        return {num:1}          #a prime number!
354           
355def makeFFTsizeList(nmin=1,nmax=1023,thresh=15):
356    ''' Provide list of optimal data sizes for FFT calculations
357
358    :param int nmin: minimum data size >= 1
359    :param int nmax: maximum data size > nmin
360    :param int thresh: maximum prime factor allowed
361    :Returns: list of data sizes where the maximum prime factor is < thresh
362    ''' 
363    plist = []
364    nmin = max(1,nmin)
365    nmax = max(nmin+1,nmax)
366    for p in range(nmin,nmax):
367        if max(factorize(p).keys()) < thresh:
368            plist.append(p)
369    return plist
370
371np.seterr(divide='ignore')
372
373# Normal distribution
374
375# loc = mu, scale = std
376_norm_pdf_C = 1./math.sqrt(2*math.pi)
377class norm_gen(st.rv_continuous):
378    'needs a doc string'
379     
380    def pdf(self,x,*args,**kwds):
381        loc,scale=kwds['loc'],kwds['scale']
382        x = (x-loc)/scale
383        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
384       
385norm = norm_gen(name='norm',longname='A normal',extradoc="""
386
387Normal distribution
388
389The location (loc) keyword specifies the mean.
390The scale (scale) keyword specifies the standard deviation.
391
392normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
393""")
394
395## Cauchy
396
397# median = loc
398
399class cauchy_gen(st.rv_continuous):
400    'needs a doc string'
401
402    def pdf(self,x,*args,**kwds):
403        loc,scale=kwds['loc'],kwds['scale']
404        x = (x-loc)/scale
405        return 1.0/np.pi/(1.0+x*x) / scale
406       
407cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
408
409Cauchy distribution
410
411cauchy.pdf(x) = 1/(pi*(1+x**2))
412
413This is the t distribution with one degree of freedom.
414""")
415   
416   
417#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
418
419
420class fcjde_gen(st.rv_continuous):
421    """
422    Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L
423    Ref: J. Appl. Cryst. (1994) 27, 892-900.
424
425    :param x: array -1 to 1
426    :param t: 2-theta position of peak
427    :param s: sum(S/L,H/L); S: sample height, H: detector opening,
428      L: sample to detector opening distance
429    :param dx: 2-theta step size in deg
430
431    :returns: for fcj.pdf
432
433     * T = x*dx+t
434     * s = S/L+H/L
435     * if x < 0::
436
437        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
438
439     * if x >= 0: fcj.pdf = 0   
440    """
441    def _pdf(self,x,t,s,dx):
442        T = dx*x+t
443        ax2 = abs(npcosd(T))
444        ax = ax2**2
445        bx = npcosd(t)**2
446        bx = np.where(ax>bx,bx,ax)
447        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
448        fx = np.where(fx > 0.,fx,0.0)
449        return fx
450             
451    def pdf(self,x,*args,**kwds):
452        loc=kwds['loc']
453        return self._pdf(x-loc,*args)
454       
455fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
456               
457def getWidthsCW(pos,sig,gam,shl):
458    '''Compute the peak widths used for computing the range of a peak
459    for constant wavelength data.
460    On low-angle side, 10 FWHM are used, on high-angle side 15 are used
461    (for peaks above 90 deg, these are reversed.)
462    '''
463    widths = [np.sqrt(sig)/100.,gam/200.]
464    fwhm = 2.355*widths[0]+2.*widths[1]
465    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
466    fmax = 15.0*fwhm
467    if pos > 90:
468        fmin,fmax = [fmax,fmin]         
469    return widths,fmin,fmax
470   
471def getWidthsTOF(pos,alp,bet,sig,gam):
472    'needs a doc string'
473    lnf = 3.3      # =log(0.001/2)
474    widths = [np.sqrt(sig),gam]
475    fwhm = 2.355*widths[0]+2.*widths[1]
476    fmin = 10.*fwhm*(1.+1./alp)   
477    fmax = 10.*fwhm*(1.+1./bet)
478    return widths,fmin,fmax
479   
480def getFWHM(pos,Inst):
481    'needs a doc string'
482    sig = lambda Th,U,V,W: 1.17741*math.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W))*math.pi/180.
483    sigTOF = lambda dsp,S0,S1,S2,Sq:  S0+S1*dsp**2+S2*dsp**4+Sq/dsp**2
484    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/180.
485    gamTOF = lambda dsp,X,Y: X*dsp+Y*dsp**2
486    if 'C' in Inst['Type'][0]:
487        s = sig(pos/2.,Inst['U'][1],Inst['V'][1],Inst['W'][1])*100.
488        g = gam(pos/2.,Inst['X'][1],Inst['Y'][1])*100.
489    else:
490        dsp = pos/Inst['difC'][0]
491        s = sigTOF(dsp,Inst['sig-0'][1],Inst['sig-1'][1],Inst['sig-2'][1],Inst['sig-q'][1])
492        g = gamTOF(dsp,Inst['X'][1],Inst['Y'][1])
493    return getgamFW(g,s)
494   
495def getgamFW(g,s):
496    'needs a doc string'
497    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.)
498    return gamFW(g,s)
499               
500def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
501    'needs a doc string'
502    DX = xdata[1]-xdata[0]
503    widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
504    x = np.linspace(pos-fmin,pos+fmin,256)
505    dx = x[1]-x[0]
506    Norm = norm.pdf(x,loc=pos,scale=widths[0])
507    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
508    arg = [pos,shl/57.2958,dx,]
509    FCJ = fcjde.pdf(x,*arg,loc=pos)
510    if len(np.nonzero(FCJ)[0])>5:
511        z = np.column_stack([Norm,Cauchy,FCJ]).T
512        Z = fft(z)
513        Df = ifft(Z.prod(axis=0)).real
514    else:
515        z = np.column_stack([Norm,Cauchy]).T
516        Z = fft(z)
517        Df = fftshift(ifft(Z.prod(axis=0))).real
518    Df /= np.sum(Df)
519    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
520    return intens*Df(xdata)*DX/dx
521
522def getBackground(pfx,parmDict,bakType,dataType,xdata):
523    'needs a doc string'
524    if 'T' in dataType:
525        q = 2.*np.pi*parmDict[pfx+'difC']/xdata
526    elif 'C' in dataType:
527        wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0))
528        q = 2.*np.pi*npsind(xdata/2.)/wave
529    yb = np.zeros_like(xdata)
530    nBak = 0
531    cw = np.diff(xdata)
532    cw = np.append(cw,cw[-1])
533    sumBk = [0.,0.,0]
534    while True:
535        key = pfx+'Back;'+str(nBak)
536        if key in parmDict:
537            nBak += 1
538        else:
539            break
540#empirical functions
541    if bakType in ['chebyschev','cosine']:
542        dt = xdata[-1]-xdata[0]   
543        for iBak in range(nBak):
544            key = pfx+'Back;'+str(iBak)
545            if bakType == 'chebyschev':
546                ybi = parmDict[key]*(2.*(xdata-xdata[0])/dt-1.)**iBak
547            elif bakType == 'cosine':
548                ybi = parmDict[key]*npcosd(xdata*iBak)
549            yb += ybi
550        sumBk[0] = np.sum(yb)
551    elif bakType in ['Q^2 power series','Q^-2 powder series']:
552        QT = 1.
553        yb += np.ones_like(yb)*parmDict[pfx+'Back;0']
554        for iBak in range(nBak-1):
555            key = pfx+'Back;'+str(iBak+1)
556            if '-2' in bakType:
557                QT *= (iBak+1)*q**-2
558            else:
559                QT *= q**2/(iBak+1)
560            yb += QT*parmDict[key]
561        sumBk[0] = np.sum(yb)
562    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
563        if nBak == 1:
564            yb = np.ones_like(xdata)*parmDict[pfx+'Back;0']
565        elif nBak == 2:
566            dX = xdata[-1]-xdata[0]
567            T2 = (xdata-xdata[0])/dX
568            T1 = 1.0-T2
569            yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2
570        else:
571            if bakType == 'lin interpolate':
572                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
573            elif bakType == 'inv interpolate':
574                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
575            elif bakType == 'log interpolate':
576                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
577            bakPos[0] = xdata[0]
578            bakPos[-1] = xdata[-1]
579            bakVals = np.zeros(nBak)
580            for i in range(nBak):
581                bakVals[i] = parmDict[pfx+'Back;'+str(i)]
582            bakInt = si.interp1d(bakPos,bakVals,'linear')
583            yb = bakInt(xdata)
584        sumBk[0] = np.sum(yb)
585#Debye function       
586    if pfx+'difC' in parmDict:
587        ff = 1.
588    else:       
589        try:
590            wave = parmDict[pfx+'Lam']
591        except KeyError:
592            wave = parmDict[pfx+'Lam1']
593        SQ = (q/(4.*np.pi))**2
594        FF = G2elem.GetFormFactorCoeff('Si')[0]
595        ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2
596    iD = 0       
597    while True:
598        try:
599            dbA = parmDict[pfx+'DebyeA;'+str(iD)]
600            dbR = parmDict[pfx+'DebyeR;'+str(iD)]
601            dbU = parmDict[pfx+'DebyeU;'+str(iD)]
602            ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
603            yb += ybi
604            sumBk[1] += np.sum(ybi)
605            iD += 1       
606        except KeyError:
607            break
608#peaks
609    iD = 0
610    while True:
611        try:
612            pkP = parmDict[pfx+'BkPkpos;'+str(iD)]
613            pkI = parmDict[pfx+'BkPkint;'+str(iD)]
614            pkS = parmDict[pfx+'BkPksig;'+str(iD)]
615            pkG = parmDict[pfx+'BkPkgam;'+str(iD)]
616            if 'C' in dataType:
617                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
618            else: #'T'OF
619                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
620            iBeg = np.searchsorted(xdata,pkP-fmin)
621            iFin = np.searchsorted(xdata,pkP+fmax)
622            lenX = len(xdata)
623            if not iBeg:
624                iFin = np.searchsorted(xdata,pkP+fmax)
625            elif iBeg == lenX:
626                iFin = iBeg
627            else:
628                iFin = np.searchsorted(xdata,pkP+fmax)
629            if 'C' in dataType:
630                ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])
631                yb[iBeg:iFin] += ybi
632            else:   #'T'OF
633                ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
634                yb[iBeg:iFin] += ybi
635            sumBk[2] += np.sum(ybi)
636            iD += 1       
637        except KeyError:
638            break
639        except ValueError:
640            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
641            break       
642    return yb,sumBk
643   
644def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata):
645    'needs a doc string'
646    if 'T' in dataType:
647        q = 2.*np.pi*parmDict[hfx+'difC']/xdata
648    elif 'C' in dataType:
649        wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0))
650        q = 2.*np.pi*npsind(xdata/2.)/wave
651    nBak = 0
652    while True:
653        key = hfx+'Back;'+str(nBak)
654        if key in parmDict:
655            nBak += 1
656        else:
657            break
658    dydb = np.zeros(shape=(nBak,len(xdata)))
659    dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata)))
660    dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata)))
661    cw = np.diff(xdata)
662    cw = np.append(cw,cw[-1])
663
664    if bakType in ['chebyschev','cosine']:
665        dt = xdata[-1]-xdata[0]   
666        for iBak in range(nBak):   
667            if bakType == 'chebyschev':
668                dydb[iBak] = (2.*(xdata-xdata[0])/dt-1.)**iBak
669            elif bakType == 'cosine':
670                dydb[iBak] = npcosd(xdata*iBak)
671    elif bakType in ['Q^2 power series','Q^-2 powder series']:
672        QT = 1.
673        dydb[0] = np.ones_like(xdata)
674        for iBak in range(nBak-1):
675            if '-2' in bakType:
676                QT *= (iBak+1)*q**-2
677            else:
678                QT *= q**2/(iBak+1)
679            dydb[iBak+1] = QT
680    elif bakType in ['lin interpolate','inv interpolate','log interpolate',]:
681        if nBak == 1:
682            dydb[0] = np.ones_like(xdata)
683        elif nBak == 2:
684            dX = xdata[-1]-xdata[0]
685            T2 = (xdata-xdata[0])/dX
686            T1 = 1.0-T2
687            dydb = [T1,T2]
688        else:
689            if bakType == 'lin interpolate':
690                bakPos = np.linspace(xdata[0],xdata[-1],nBak,True)
691            elif bakType == 'inv interpolate':
692                bakPos = 1./np.linspace(1./xdata[-1],1./xdata[0],nBak,True)
693            elif bakType == 'log interpolate':
694                bakPos = np.exp(np.linspace(np.log(xdata[0]),np.log(xdata[-1]),nBak,True))
695            bakPos[0] = xdata[0]
696            bakPos[-1] = xdata[-1]
697            for i,pos in enumerate(bakPos):
698                if i == 0:
699                    dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.)
700                elif i == len(bakPos)-1:
701                    dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.)
702                else:
703                    dydb[i] = np.where(xdata>bakPos[i],
704                        np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.),
705                        np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.))
706    if hfx+'difC' in parmDict:
707        ff = 1.
708    else:
709        try:
710            wave = parmDict[hfx+'Lam']
711        except KeyError:
712            wave = parmDict[hfx+'Lam1']
713        q = 4.0*np.pi*npsind(xdata/2.0)/wave
714        SQ = (q/(4*np.pi))**2
715        FF = G2elem.GetFormFactorCoeff('Si')[0]
716        ff = np.array(G2elem.ScatFac(FF,SQ)[0])
717    iD = 0       
718    while True:
719        try:
720            if hfx+'difC' in parmDict:
721                q = 2*np.pi*parmDict[hfx+'difC']/xdata
722            dbA = parmDict[hfx+'DebyeA;'+str(iD)]
723            dbR = parmDict[hfx+'DebyeR;'+str(iD)]
724            dbU = parmDict[hfx+'DebyeU;'+str(iD)]
725            sqr = np.sin(q*dbR)/(q*dbR)
726            cqr = np.cos(q*dbR)
727            temp = np.exp(-dbU*q**2)
728            dyddb[3*iD] = ff*sqr*temp/(np.pi*cw)
729            dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(np.pi*dbR*cw)
730            dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2/(np.pi*cw)
731            iD += 1       #ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR)
732        except KeyError:
733            break
734    iD = 0
735    while True:
736        try:
737            pkP = parmDict[hfx+'BkPkpos;'+str(iD)]
738            pkI = parmDict[hfx+'BkPkint;'+str(iD)]
739            pkS = parmDict[hfx+'BkPksig;'+str(iD)]
740            pkG = parmDict[hfx+'BkPkgam;'+str(iD)]
741            if 'C' in dataType:
742                Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002)
743            else: #'T'OF
744                Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG)
745            iBeg = np.searchsorted(xdata,pkP-fmin)
746            iFin = np.searchsorted(xdata,pkP+fmax)
747            lenX = len(xdata)
748            if not iBeg:
749                iFin = np.searchsorted(xdata,pkP+fmax)
750            elif iBeg == lenX:
751                iFin = iBeg
752            else:
753                iFin = np.searchsorted(xdata,pkP+fmax)
754            if 'C' in dataType:
755                Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin])
756                dydpk[4*iD][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdp
757                dydpk[4*iD+1][iBeg:iFin] += 100.*cw[iBeg:iFin]*Df
758                dydpk[4*iD+2][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFds
759                dydpk[4*iD+3][iBeg:iFin] += 100.*cw[iBeg:iFin]*pkI*dFdg
760            else:   #'T'OF
761                Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])
762                dydpk[4*iD][iBeg:iFin] += pkI*dFdp
763                dydpk[4*iD+1][iBeg:iFin] += Df
764                dydpk[4*iD+2][iBeg:iFin] += pkI*dFds
765                dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg
766            iD += 1       
767        except KeyError:
768            break
769        except ValueError:
770            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
771            break       
772    return dydb,dyddb,dydpk
773
774#use old fortran routine
775def getFCJVoigt3(pos,sig,gam,shl,xdata):
776    'needs a doc string'
777   
778    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
779#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
780    Df /= np.sum(Df)
781    return Df
782
783def getdFCJVoigt3(pos,sig,gam,shl,xdata):
784    'needs a doc string'
785   
786    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
787#    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
788    sumDf = np.sum(Df)
789    return Df,dFdp,dFds,dFdg,dFdsh
790
791def getPsVoigt(pos,sig,gam,xdata):
792    'needs a doc string'
793   
794    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
795    Df /= np.sum(Df)
796    return Df
797
798def getdPsVoigt(pos,sig,gam,xdata):
799    'needs a doc string'
800   
801    Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam)
802    sumDf = np.sum(Df)
803    return Df,dFdp,dFds,dFdg
804
805def getEpsVoigt(pos,alp,bet,sig,gam,xdata):
806    'needs a doc string'
807    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
808    Df /= np.sum(Df)
809    return Df 
810   
811def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
812    'needs a doc string'
813    Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
814    sumDf = np.sum(Df)
815    return Df,dFdp,dFda,dFdb,dFds,dFdg   
816
817def ellipseSize(H,Sij,GB):
818    'needs a doc string'
819    HX = np.inner(H.T,GB)
820    lenHX = np.sqrt(np.sum(HX**2))
821    Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij))           
822    R = np.inner(HX/lenHX,Rsize)*Esize         #want column length for hkl in crystal
823    lenR = np.sqrt(np.sum(R**2))
824    return lenR
825
826def ellipseSizeDerv(H,Sij,GB):
827    'needs a doc string'
828    lenR = ellipseSize(H,Sij,GB)
829    delt = 0.001
830    dRdS = np.zeros(6)
831    for i in range(6):
832        Sij[i] -= delt
833        lenM = ellipseSize(H,Sij,GB)
834        Sij[i] += 2.*delt
835        lenP = ellipseSize(H,Sij,GB)
836        Sij[i] -= delt
837        dRdS[i] = (lenP-lenM)/(2.*delt)
838    return lenR,dRdS
839
840def getHKLpeak(dmin,SGData,A,Inst=None):
841    'needs a doc string'
842    HKL = G2lat.GenHLaue(dmin,SGData,A)       
843    HKLs = []
844    for h,k,l,d in HKL:
845        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
846        if not ext:
847            if Inst == None:
848                HKLs.append([h,k,l,d,0,-1])
849            else:
850                HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1])
851    return HKLs
852
853def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A):
854    'needs a doc string'
855    HKLs = []
856    vec = np.array(Vec)
857    vstar = np.sqrt(G2lat.calc_rDsq(vec,A))     #find extra needed for -n SS reflections
858    dvec = 1./(maxH*vstar+1./dmin)
859    HKL = G2lat.GenHLaue(dvec,SGData,A)       
860    SSdH = [vec*h for h in range(-maxH,maxH+1)]
861    SSdH = dict(zip(range(-maxH,maxH+1),SSdH))
862    for h,k,l,d in HKL:
863        ext = G2spc.GenHKLf([h,k,l],SGData)[0]
864        if not ext and d >= dmin:
865            HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1])
866        for dH in SSdH:
867            if dH:
868                DH = SSdH[dH]
869                H = [h+DH[0],k+DH[1],l+DH[2]]
870                d = 1/np.sqrt(G2lat.calc_rDsq(H,A))
871                if d >= dmin:
872                    HKLM = np.array([h,k,l,dH])
873                    if G2spc.checkSSextc(HKLM,SSGData):
874                        HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1])   
875    return G2lat.sortHKLd(HKLs,True,True,True)
876
877def getPeakProfile(dataType,parmDict,xdata,varyList,bakType):
878    'needs a doc string'
879   
880    yb = getBackground('',parmDict,bakType,dataType,xdata)[0]
881    yc = np.zeros_like(yb)
882    cw = np.diff(xdata)
883    cw = np.append(cw,cw[-1])
884    if 'C' in dataType:
885        shl = max(parmDict['SH/L'],0.002)
886        Ka2 = False
887        if 'Lam1' in parmDict.keys():
888            Ka2 = True
889            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
890            kRatio = parmDict['I(L2)/I(L1)']
891        iPeak = 0
892        while True:
893            try:
894                pos = parmDict['pos'+str(iPeak)]
895                tth = (pos-parmDict['Zero'])
896                intens = parmDict['int'+str(iPeak)]
897                sigName = 'sig'+str(iPeak)
898                if sigName in varyList:
899                    sig = parmDict[sigName]
900                else:
901                    sig = G2mth.getCWsig(parmDict,tth)
902                sig = max(sig,0.001)          #avoid neg sigma
903                gamName = 'gam'+str(iPeak)
904                if gamName in varyList:
905                    gam = parmDict[gamName]
906                else:
907                    gam = G2mth.getCWgam(parmDict,tth)
908                gam = max(gam,0.001)             #avoid neg gamma
909                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
910                iBeg = np.searchsorted(xdata,pos-fmin)
911                iFin = np.searchsorted(xdata,pos+fmin)
912                if not iBeg+iFin:       #peak below low limit
913                    iPeak += 1
914                    continue
915                elif not iBeg-iFin:     #peak above high limit
916                    return yb+yc
917                yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
918                if Ka2:
919                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
920                    iBeg = np.searchsorted(xdata,pos2-fmin)
921                    iFin = np.searchsorted(xdata,pos2+fmin)
922                    if iBeg-iFin:
923                        yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
924                iPeak += 1
925            except KeyError:        #no more peaks to process
926                return yb+yc
927    else:
928        Pdabc = parmDict['Pdabc']
929        difC = parmDict['difC']
930        iPeak = 0
931        while True:
932            try:
933                pos = parmDict['pos'+str(iPeak)]               
934                tof = pos-parmDict['Zero']
935                dsp = tof/difC
936                intens = parmDict['int'+str(iPeak)]
937                alpName = 'alp'+str(iPeak)
938                if alpName in varyList:
939                    alp = parmDict[alpName]
940                else:
941                    if len(Pdabc):
942                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
943                    else:
944                        alp = G2mth.getTOFalpha(parmDict,dsp)
945                alp = max(0.0001,alp)
946                betName = 'bet'+str(iPeak)
947                if betName in varyList:
948                    bet = parmDict[betName]
949                else:
950                    if len(Pdabc):
951                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
952                    else:
953                        bet = G2mth.getTOFbeta(parmDict,dsp)
954                bet = max(0.0001,bet)
955                sigName = 'sig'+str(iPeak)
956                if sigName in varyList:
957                    sig = parmDict[sigName]
958                else:
959                    sig = G2mth.getTOFsig(parmDict,dsp)
960                gamName = 'gam'+str(iPeak)
961                if gamName in varyList:
962                    gam = parmDict[gamName]
963                else:
964                    gam = G2mth.getTOFgamma(parmDict,dsp)
965                gam = max(gam,0.001)             #avoid neg gamma
966                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
967                iBeg = np.searchsorted(xdata,pos-fmin)
968                iFin = np.searchsorted(xdata,pos+fmax)
969                lenX = len(xdata)
970                if not iBeg:
971                    iFin = np.searchsorted(xdata,pos+fmax)
972                elif iBeg == lenX:
973                    iFin = iBeg
974                else:
975                    iFin = np.searchsorted(xdata,pos+fmax)
976                if not iBeg+iFin:       #peak below low limit
977                    iPeak += 1
978                    continue
979                elif not iBeg-iFin:     #peak above high limit
980                    return yb+yc
981                yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
982                iPeak += 1
983            except KeyError:        #no more peaks to process
984                return yb+yc
985           
986def getPeakProfileDerv(dataType,parmDict,xdata,varyList,bakType):
987    'needs a doc string'
988# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
989    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
990    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,dataType,xdata)
991    if 'Back;0' in varyList:            #background derivs are in front if present
992        dMdv[0:len(dMdb)] = dMdb
993    names = ['DebyeA','DebyeR','DebyeU']
994    for name in varyList:
995        if 'Debye' in name:
996            parm,id = name.split(';')
997            ip = names.index(parm)
998            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
999    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1000    for name in varyList:
1001        if 'BkPk' in name:
1002            parm,id = name.split(';')
1003            ip = names.index(parm)
1004            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
1005    cw = np.diff(xdata)
1006    cw = np.append(cw,cw[-1])
1007    if 'C' in dataType:
1008        shl = max(parmDict['SH/L'],0.002)
1009        Ka2 = False
1010        if 'Lam1' in parmDict.keys():
1011            Ka2 = True
1012            lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
1013            kRatio = parmDict['I(L2)/I(L1)']
1014        iPeak = 0
1015        while True:
1016            try:
1017                pos = parmDict['pos'+str(iPeak)]
1018                tth = (pos-parmDict['Zero'])
1019                intens = parmDict['int'+str(iPeak)]
1020                sigName = 'sig'+str(iPeak)
1021                if sigName in varyList:
1022                    sig = parmDict[sigName]
1023                    dsdU = dsdV = dsdW = 0
1024                else:
1025                    sig = G2mth.getCWsig(parmDict,tth)
1026                    dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth)
1027                sig = max(sig,0.001)          #avoid neg sigma
1028                gamName = 'gam'+str(iPeak)
1029                if gamName in varyList:
1030                    gam = parmDict[gamName]
1031                    dgdX = dgdY = 0
1032                else:
1033                    gam = G2mth.getCWgam(parmDict,tth)
1034                    dgdX,dgdY = G2mth.getCWgamDeriv(tth)
1035                gam = max(gam,0.001)             #avoid neg gamma
1036                Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl)
1037                iBeg = np.searchsorted(xdata,pos-fmin)
1038                iFin = np.searchsorted(xdata,pos+fmin)
1039                if not iBeg+iFin:       #peak below low limit
1040                    iPeak += 1
1041                    continue
1042                elif not iBeg-iFin:     #peak above high limit
1043                    break
1044                dMdpk = np.zeros(shape=(6,len(xdata)))
1045                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
1046                for i in range(1,5):
1047                    dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]
1048                dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]
1049                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1050                if Ka2:
1051                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
1052                    iBeg = np.searchsorted(xdata,pos2-fmin)
1053                    iFin = np.searchsorted(xdata,pos2+fmin)
1054                    if iBeg-iFin:
1055                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
1056                        for i in range(1,5):
1057                            dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]
1058                        dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]
1059                        dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]
1060                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
1061                for parmName in ['pos','int','sig','gam']:
1062                    try:
1063                        idx = varyList.index(parmName+str(iPeak))
1064                        dMdv[idx] = dervDict[parmName]
1065                    except ValueError:
1066                        pass
1067                if 'U' in varyList:
1068                    dMdv[varyList.index('U')] += dsdU*dervDict['sig']
1069                if 'V' in varyList:
1070                    dMdv[varyList.index('V')] += dsdV*dervDict['sig']
1071                if 'W' in varyList:
1072                    dMdv[varyList.index('W')] += dsdW*dervDict['sig']
1073                if 'X' in varyList:
1074                    dMdv[varyList.index('X')] += dgdX*dervDict['gam']
1075                if 'Y' in varyList:
1076                    dMdv[varyList.index('Y')] += dgdY*dervDict['gam']
1077                if 'SH/L' in varyList:
1078                    dMdv[varyList.index('SH/L')] += dervDict['shl']         #problem here
1079                if 'I(L2)/I(L1)' in varyList:
1080                    dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2']
1081                iPeak += 1
1082            except KeyError:        #no more peaks to process
1083                break
1084    else:
1085        Pdabc = parmDict['Pdabc']
1086        difC = parmDict['difC']
1087        iPeak = 0
1088        while True:
1089            try:
1090                pos = parmDict['pos'+str(iPeak)]               
1091                tof = pos-parmDict['Zero']
1092                dsp = tof/difC
1093                intens = parmDict['int'+str(iPeak)]
1094                alpName = 'alp'+str(iPeak)
1095                if alpName in varyList:
1096                    alp = parmDict[alpName]
1097                else:
1098                    if len(Pdabc):
1099                        alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1100                        dad0 = 0
1101                    else:
1102                        alp = G2mth.getTOFalpha(parmDict,dsp)
1103                        dada0 = G2mth.getTOFalphaDeriv(dsp)
1104                betName = 'bet'+str(iPeak)
1105                if betName in varyList:
1106                    bet = parmDict[betName]
1107                else:
1108                    if len(Pdabc):
1109                        bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1110                        dbdb0 = dbdb1 = dbdb2 = 0
1111                    else:
1112                        bet = G2mth.getTOFbeta(parmDict,dsp)
1113                        dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp)
1114                sigName = 'sig'+str(iPeak)
1115                if sigName in varyList:
1116                    sig = parmDict[sigName]
1117                    dsds0 = dsds1 = dsds2 = dsds3 = 0
1118                else:
1119                    sig = G2mth.getTOFsig(parmDict,dsp)
1120                    dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp)
1121                gamName = 'gam'+str(iPeak)
1122                if gamName in varyList:
1123                    gam = parmDict[gamName]
1124                    dsdX = dsdY = 0
1125                else:
1126                    gam = G2mth.getTOFgamma(parmDict,dsp)
1127                    dsdX,dsdY = G2mth.getTOFgammaDeriv(dsp)
1128                gam = max(gam,0.001)             #avoid neg gamma
1129                Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam)
1130                iBeg = np.searchsorted(xdata,pos-fmin)
1131                lenX = len(xdata)
1132                if not iBeg:
1133                    iFin = np.searchsorted(xdata,pos+fmax)
1134                elif iBeg == lenX:
1135                    iFin = iBeg
1136                else:
1137                    iFin = np.searchsorted(xdata,pos+fmax)
1138                if not iBeg+iFin:       #peak below low limit
1139                    iPeak += 1
1140                    continue
1141                elif not iBeg-iFin:     #peak above high limit
1142                    break
1143                dMdpk = np.zeros(shape=(7,len(xdata)))
1144                dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])
1145                for i in range(1,6):
1146                    dMdpk[i][iBeg:iFin] += intens*cw[iBeg:iFin]*dMdipk[i]
1147                dMdpk[0][iBeg:iFin] += cw[iBeg:iFin]*dMdipk[0]
1148                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}
1149                for parmName in ['pos','int','alp','bet','sig','gam']:
1150                    try:
1151                        idx = varyList.index(parmName+str(iPeak))
1152                        dMdv[idx] = dervDict[parmName]
1153                    except ValueError:
1154                        pass
1155                if 'alpha' in varyList:
1156                    dMdv[varyList.index('alpha')] += dada0*dervDict['alp']
1157                if 'beta-0' in varyList:
1158                    dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet']
1159                if 'beta-1' in varyList:
1160                    dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet']
1161                if 'beta-q' in varyList:
1162                    dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet']
1163                if 'sig-0' in varyList:
1164                    dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig']
1165                if 'sig-1' in varyList:
1166                    dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig']
1167                if 'sig-2' in varyList:
1168                    dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig']
1169                if 'sig-q' in varyList:
1170                    dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig']
1171                if 'X' in varyList:
1172                    dMdv[varyList.index('X')] += dsdX*dervDict['gam']
1173                if 'Y' in varyList:
1174                    dMdv[varyList.index('Y')] += dsdY*dervDict['gam']         #problem here
1175                iPeak += 1
1176            except KeyError:        #no more peaks to process
1177                break
1178    return dMdv
1179       
1180def Dict2Values(parmdict, varylist):
1181    '''Use before call to leastsq to setup list of values for the parameters
1182    in parmdict, as selected by key in varylist'''
1183    return [parmdict[key] for key in varylist] 
1184   
1185def Values2Dict(parmdict, varylist, values):
1186    ''' Use after call to leastsq to update the parameter dictionary with
1187    values corresponding to keys in varylist'''
1188    parmdict.update(zip(varylist,values))
1189   
1190def SetBackgroundParms(Background):
1191    'needs a doc string'
1192    if len(Background) == 1:            # fix up old backgrounds
1193        BackGround.append({'nDebye':0,'debyeTerms':[]})
1194    bakType,bakFlag = Background[0][:2]
1195    backVals = Background[0][3:]
1196    backNames = ['Back;'+str(i) for i in range(len(backVals))]
1197    Debye = Background[1]           #also has background peaks stuff
1198    backDict = dict(zip(backNames,backVals))
1199    backVary = []
1200    if bakFlag:
1201        backVary = backNames
1202
1203    backDict['nDebye'] = Debye['nDebye']
1204    debyeDict = {}
1205    debyeList = []
1206    for i in range(Debye['nDebye']):
1207        debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)]
1208        debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1209        debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1210    debyeVary = []
1211    for item in debyeList:
1212        if item[1]:
1213            debyeVary.append(item[0])
1214    backDict.update(debyeDict)
1215    backVary += debyeVary
1216
1217    backDict['nPeaks'] = Debye['nPeaks']
1218    peaksDict = {}
1219    peaksList = []
1220    for i in range(Debye['nPeaks']):
1221        peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)]
1222        peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2])))
1223        peaksList += zip(peaksNames,Debye['peaksList'][i][1::2])
1224    peaksVary = []
1225    for item in peaksList:
1226        if item[1]:
1227            peaksVary.append(item[0])
1228    backDict.update(peaksDict)
1229    backVary += peaksVary   
1230    return bakType,backDict,backVary
1231   
1232def DoCalibInst(IndexPeaks,Inst):
1233   
1234    def SetInstParms():
1235        dataType = Inst['Type'][0]
1236        insVary = []
1237        insNames = []
1238        insVals = []
1239        for parm in Inst:
1240            insNames.append(parm)
1241            insVals.append(Inst[parm][1])
1242            if parm in ['Lam','difC','difA','difB','Zero',]:
1243                if Inst[parm][2]:
1244                    insVary.append(parm)
1245        instDict = dict(zip(insNames,insVals))
1246        return dataType,instDict,insVary
1247       
1248    def GetInstParms(parmDict,Inst,varyList):
1249        for name in Inst:
1250            Inst[name][1] = parmDict[name]
1251       
1252    def InstPrint(Inst,sigDict):
1253        print 'Instrument Parameters:'
1254        if 'C' in Inst['Type'][0]:
1255            ptfmt = "%12.6f"
1256        else:
1257            ptfmt = "%12.3f"
1258        ptlbls = 'names :'
1259        ptstr =  'values:'
1260        sigstr = 'esds  :'
1261        for parm in Inst:
1262            if parm in  ['Lam','difC','difA','difB','Zero',]:
1263                ptlbls += "%s" % (parm.center(12))
1264                ptstr += ptfmt % (Inst[parm][1])
1265                if parm in sigDict:
1266                    sigstr += ptfmt % (sigDict[parm])
1267                else:
1268                    sigstr += 12*' '
1269        print ptlbls
1270        print ptstr
1271        print sigstr
1272       
1273    def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList):
1274        parmDict.update(zip(varyList,values))
1275        return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos)
1276
1277    peakPos = []
1278    peakDsp = []
1279    peakWt = []
1280    for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]):
1281        if peak[2] and peak[3]:
1282            peakPos.append(peak[0])
1283            peakDsp.append(peak[-1])    #d-calc
1284            peakWt.append(peak[-1]**2/sig**2)   #weight by d**2
1285    peakPos = np.array(peakPos)
1286    peakDsp = np.array(peakDsp)
1287    peakWt = np.array(peakWt)
1288    dataType,insDict,insVary = SetInstParms()
1289    parmDict = {}
1290    parmDict.update(insDict)
1291    varyList = insVary
1292    if not len(varyList):
1293        print '**** ERROR - nothing to refine! ****'
1294        return False
1295    while True:
1296        begin = time.time()
1297        values =  np.array(Dict2Values(parmDict, varyList))
1298        result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001,
1299            args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList))
1300        ncyc = int(result[2]['nfev']/2)
1301        runtime = time.time()-begin   
1302        chisq = np.sum(result[2]['fvec']**2)
1303        Values2Dict(parmDict, varyList, result[0])
1304        GOF = chisq/(len(peakPos)-len(varyList))       #reduced chi^2
1305        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(peakPos),' Number of parameters: ',len(varyList)
1306        print 'calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1307        print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1308        try:
1309            sig = np.sqrt(np.diag(result[1])*GOF)
1310            if np.any(np.isnan(sig)):
1311                print '*** Least squares aborted - some invalid esds possible ***'
1312            break                   #refinement succeeded - finish up!
1313        except ValueError:          #result[1] is None on singular matrix
1314            print '**** Refinement failed - singular matrix ****'
1315       
1316    sigDict = dict(zip(varyList,sig))
1317    GetInstParms(parmDict,Inst,varyList)
1318    InstPrint(Inst,sigDict)
1319    return True
1320           
1321def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,prevVaryList=[],oneCycle=False,controls=None,dlg=None):
1322    'needs a doc string'
1323       
1324    def GetBackgroundParms(parmList,Background):
1325        iBak = 0
1326        while True:
1327            try:
1328                bakName = 'Back;'+str(iBak)
1329                Background[0][iBak+3] = parmList[bakName]
1330                iBak += 1
1331            except KeyError:
1332                break
1333        iDb = 0
1334        while True:
1335            names = ['DebyeA;','DebyeR;','DebyeU;']
1336            try:
1337                for i,name in enumerate(names):
1338                    val = parmList[name+str(iDb)]
1339                    Background[1]['debyeTerms'][iDb][2*i] = val
1340                iDb += 1
1341            except KeyError:
1342                break
1343        iDb = 0
1344        while True:
1345            names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;']
1346            try:
1347                for i,name in enumerate(names):
1348                    val = parmList[name+str(iDb)]
1349                    Background[1]['peaksList'][iDb][2*i] = val
1350                iDb += 1
1351            except KeyError:
1352                break
1353               
1354    def BackgroundPrint(Background,sigDict):
1355        if Background[0][1]:
1356            print 'Background coefficients for',Background[0][0],'function'
1357            ptfmt = "%12.5f"
1358            ptstr =  'values:'
1359            sigstr = 'esds  :'
1360            for i,back in enumerate(Background[0][3:]):
1361                ptstr += ptfmt % (back)
1362                sigstr += ptfmt % (sigDict['Back;'+str(i)])
1363            print ptstr
1364            print sigstr
1365        else:
1366            print 'Background not refined'
1367        if Background[1]['nDebye']:
1368            parms = ['DebyeA;','DebyeR;','DebyeU;']
1369            print 'Debye diffuse scattering coefficients'
1370            ptfmt = "%12.5f"
1371            print ' term       DebyeA       esd        DebyeR       esd        DebyeU        esd'
1372            for term in range(Background[1]['nDebye']):
1373                line = ' term %d'%(term)
1374                for ip,name in enumerate(parms):
1375                    line += ptfmt%(Background[1]['debyeTerms'][term][2*ip])
1376                    if name+str(term) in sigDict:
1377                        line += ptfmt%(sigDict[name+str(term)])
1378                print line
1379        if Background[1]['nPeaks']:
1380            parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
1381            print 'Peaks in background coefficients'
1382            ptfmt = "%15.3f"
1383            names =   'names :'
1384            ptstr =  'values:'
1385            sigstr = 'esds  :'
1386            for item in sigDict:
1387                if 'BkPk' in item:
1388                    names += '%15s'%(item)
1389                    sigstr += ptfmt%(sigDict[item])
1390                    parm,id = item.split(';')
1391                    ip = parms.index(parm)
1392                    ptstr += ptfmt%(Background[1]['peaksList'][int(id)][2*ip])
1393            print names
1394            print ptstr
1395            print sigstr
1396                           
1397    def SetInstParms(Inst):
1398        dataType = Inst['Type'][0]
1399        insVary = []
1400        insNames = []
1401        insVals = []
1402        for parm in Inst:
1403            insNames.append(parm)
1404            insVals.append(Inst[parm][1])
1405            if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1406                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1407                    insVary.append(parm)
1408        instDict = dict(zip(insNames,insVals))
1409        instDict['X'] = max(instDict['X'],0.01)
1410        instDict['Y'] = max(instDict['Y'],0.01)
1411        if 'SH/L' in instDict:
1412            instDict['SH/L'] = max(instDict['SH/L'],0.002)
1413        return dataType,instDict,insVary
1414       
1415    def GetInstParms(parmDict,Inst,varyList,Peaks):
1416        for name in Inst:
1417            Inst[name][1] = parmDict[name]
1418        iPeak = 0
1419        while True:
1420            try:
1421                sigName = 'sig'+str(iPeak)
1422                pos = parmDict['pos'+str(iPeak)]
1423                if sigName not in varyList:
1424                    if 'C' in Inst['Type'][0]:
1425                        parmDict[sigName] = G2mth.getCWsig(parmDict,pos)
1426                    else:
1427                        dsp = G2lat.Pos2dsp(Inst,pos)
1428                        parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp)
1429                gamName = 'gam'+str(iPeak)
1430                if gamName not in varyList:
1431                    if 'C' in Inst['Type'][0]:
1432                        parmDict[gamName] = G2mth.getCWgam(parmDict,pos)
1433                    else:
1434                        dsp = G2lat.Pos2dsp(Inst,pos)
1435                        parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp)
1436                iPeak += 1
1437            except KeyError:
1438                break
1439       
1440    def InstPrint(Inst,sigDict):
1441        print 'Instrument Parameters:'
1442        ptfmt = "%12.6f"
1443        ptlbls = 'names :'
1444        ptstr =  'values:'
1445        sigstr = 'esds  :'
1446        for parm in Inst:
1447            if parm in  ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1448                'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',]:
1449                ptlbls += "%s" % (parm.center(12))
1450                ptstr += ptfmt % (Inst[parm][1])
1451                if parm in sigDict:
1452                    sigstr += ptfmt % (sigDict[parm])
1453                else:
1454                    sigstr += 12*' '
1455        print ptlbls
1456        print ptstr
1457        print sigstr
1458
1459    def SetPeaksParms(dataType,Peaks):
1460        peakNames = []
1461        peakVary = []
1462        peakVals = []
1463        if 'C' in dataType:
1464            names = ['pos','int','sig','gam']
1465        else:
1466            names = ['pos','int','alp','bet','sig','gam']
1467        for i,peak in enumerate(Peaks):
1468            for j,name in enumerate(names):
1469                peakVals.append(peak[2*j])
1470                parName = name+str(i)
1471                peakNames.append(parName)
1472                if peak[2*j+1]:
1473                    peakVary.append(parName)
1474        return dict(zip(peakNames,peakVals)),peakVary
1475               
1476    def GetPeaksParms(Inst,parmDict,Peaks,varyList):
1477        if 'C' in Inst['Type'][0]:
1478            names = ['pos','int','sig','gam']
1479        else:   #'T'
1480            names = ['pos','int','alp','bet','sig','gam']
1481        for i,peak in enumerate(Peaks):
1482            pos = parmDict['pos'+str(i)]
1483            if 'difC' in Inst:
1484                dsp = pos/Inst['difC'][1]
1485            for j in range(len(names)):
1486                parName = names[j]+str(i)
1487                if parName in varyList:
1488                    peak[2*j] = parmDict[parName]
1489                elif 'alpha' in parName:
1490                    peak[2*j] = parmDict['alpha']/dsp
1491                elif 'beta' in parName:
1492                    peak[2*j] = G2mth.getTOFbeta(parmDict,dsp)
1493                elif 'sig' in parName:
1494                    if 'C' in Inst['Type'][0]:
1495                        peak[2*j] = G2mth.getCWsig(parmDict,pos)
1496                    else:
1497                        peak[2*j] = G2mth.getTOFsig(parmDict,dsp)
1498                elif 'gam' in parName:
1499                    if 'C' in Inst['Type'][0]:
1500                        peak[2*j] = G2mth.getCWgam(parmDict,pos)
1501                    else:
1502                        peak[2*j] = G2mth.getTOFgamma(parmDict,dsp)
1503                       
1504    def PeaksPrint(dataType,parmDict,sigDict,varyList):
1505        print 'Peak coefficients:'
1506        if 'C' in dataType:
1507            names = ['pos','int','sig','gam']
1508        else:   #'T'
1509            names = ['pos','int','alp','bet','sig','gam']           
1510        head = 13*' '
1511        for name in names:
1512            if name in ['alp','bet']:
1513                head += name.center(8)+'esd'.center(8)
1514            else:
1515                head += name.center(10)+'esd'.center(10)
1516        print head
1517        if 'C' in dataType:
1518            ptfmt = {'pos':"%10.5f",'int':"%10.1f",'sig':"%10.3f",'gam':"%10.3f"}
1519        else:
1520            ptfmt = {'pos':"%10.2f",'int':"%10.4f",'alp':"%8.3f",'bet':"%8.5f",'sig':"%10.3f",'gam':"%10.3f"}
1521        for i,peak in enumerate(Peaks):
1522            ptstr =  ':'
1523            for j in range(len(names)):
1524                name = names[j]
1525                parName = name+str(i)
1526                ptstr += ptfmt[name] % (parmDict[parName])
1527                if parName in varyList:
1528                    ptstr += ptfmt[name] % (sigDict[parName])
1529                else:
1530                    if name in ['alp','bet']:
1531                        ptstr += 8*' '
1532                    else:
1533                        ptstr += 10*' '
1534            print '%s'%(('Peak'+str(i+1)).center(8)),ptstr
1535               
1536    def devPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):
1537        parmdict.update(zip(varylist,values))
1538        return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,varylist,bakType)
1539           
1540    def errPeakProfile(values,xdata,ydata, weights,dataType,parmdict,varylist,bakType,dlg):       
1541        parmdict.update(zip(varylist,values))
1542        M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,varylist,bakType)-ydata)
1543        Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.)
1544        if dlg:
1545            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0]
1546            if not GoOn:
1547                return -M           #abort!!
1548        return M
1549       
1550    if controls:
1551        Ftol = controls['min dM/M']
1552        derivType = controls['deriv type']
1553    else:
1554        Ftol = 0.0001
1555        derivType = 'analytic'
1556    if oneCycle:
1557        Ftol = 1.0
1558    x,y,w,yc,yb,yd = data               #these are numpy arrays!
1559    yc *= 0.                            #set calcd ones to zero
1560    yb *= 0.
1561    yd *= 0.
1562    xBeg = np.searchsorted(x,Limits[0])
1563    xFin = np.searchsorted(x,Limits[1])
1564    bakType,bakDict,bakVary = SetBackgroundParms(Background)
1565    dataType,insDict,insVary = SetInstParms(Inst)
1566    peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks)
1567    parmDict = {}
1568    parmDict.update(bakDict)
1569    parmDict.update(insDict)
1570    parmDict.update(peakDict)
1571    parmDict['Pdabc'] = []      #dummy Pdabc
1572    parmDict.update(Inst2)      #put in real one if there
1573    if prevVaryList:
1574        varyList = prevVaryList[:]
1575    else:
1576        varyList = bakVary+insVary+peakVary
1577    fullvaryList = varyList[:]
1578    while True:
1579        begin = time.time()
1580        values =  np.array(Dict2Values(parmDict, varyList))
1581        Rvals = {}
1582        badVary = []
1583        if FitPgm == 'LSQ':
1584            try:
1585                result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True,
1586                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg))
1587                ncyc = int(result[2]['nfev']/2)
1588            finally:
1589                dlg.Destroy()
1590            runtime = time.time()-begin   
1591            chisq = np.sum(result[2]['fvec']**2)
1592            Values2Dict(parmDict, varyList, result[0])
1593            Rvals['Rwp'] = np.sqrt(chisq/np.sum(w[xBeg:xFin]*y[xBeg:xFin]**2))*100.      #to %
1594            Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList))       #reduced chi^2
1595            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',xFin-xBeg,' Number of parameters: ',len(varyList)
1596            print 'fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1597            print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1598            try:
1599                sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
1600                if np.any(np.isnan(sig)):
1601                    print '*** Least squares aborted - some invalid esds possible ***'
1602                break                   #refinement succeeded - finish up!
1603            except ValueError:          #result[1] is None on singular matrix
1604                print '**** Refinement failed - singular matrix ****'
1605                Ipvt = result[2]['ipvt']
1606                for i,ipvt in enumerate(Ipvt):
1607                    if not np.sum(result[2]['fjac'],axis=1)[i]:
1608                        print 'Removing parameter: ',varyList[ipvt-1]
1609                        badVary.append(varyList[ipvt-1])
1610                        del(varyList[ipvt-1])
1611                        break
1612        elif FitPgm == 'BFGS':
1613            print 'Other program here'
1614            return {}
1615       
1616    sigDict = dict(zip(varyList,sig))
1617    yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin])[0]
1618    yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],varyList,bakType)
1619    yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin]
1620    GetBackgroundParms(parmDict,Background)
1621    BackgroundPrint(Background,sigDict)
1622    GetInstParms(parmDict,Inst,varyList,Peaks)
1623    InstPrint(Inst,sigDict)
1624    GetPeaksParms(Inst,parmDict,Peaks,varyList)   
1625    PeaksPrint(dataType,parmDict,sigDict,varyList)
1626    return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary
1627
1628def calcIncident(Iparm,xdata):
1629    'needs a doc string'
1630
1631    def IfunAdv(Iparm,xdata):
1632        Itype = Iparm['Itype']
1633        Icoef = Iparm['Icoeff']
1634        DYI = np.ones((12,xdata.shape[0]))
1635        YI = np.ones_like(xdata)*Icoef[0]
1636       
1637        x = xdata/1000.                 #expressions are in ms
1638        if Itype == 'Exponential':
1639            for i in [1,3,5,7,9]:
1640                Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1641                YI += Icoef[i]*Eterm
1642                DYI[i] *= Eterm
1643                DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)           
1644        elif 'Maxwell'in Itype:
1645            Eterm = np.exp(-Icoef[2]/x**2)
1646            DYI[1] = Eterm/x**5
1647            DYI[2] = -Icoef[1]*DYI[1]/x**2
1648            YI += (Icoef[1]*Eterm/x**5)
1649            if 'Exponential' in Itype:
1650                for i in range(3,12,2):
1651                    Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2))
1652                    YI += Icoef[i]*Eterm
1653                    DYI[i] *= Eterm
1654                    DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2)
1655            else:   #Chebyschev
1656                T = (2./x)-1.
1657                Ccof = np.ones((12,xdata.shape[0]))
1658                Ccof[1] = T
1659                for i in range(2,12):
1660                    Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2]
1661                for i in range(1,10):
1662                    YI += Ccof[i]*Icoef[i+2]
1663                    DYI[i+2] =Ccof[i]
1664        return YI,DYI
1665       
1666    Iesd = np.array(Iparm['Iesd'])
1667    Icovar = Iparm['Icovar']
1668    YI,DYI = IfunAdv(Iparm,xdata)
1669    YI = np.where(YI>0,YI,1.)
1670    WYI = np.zeros_like(xdata)
1671    vcov = np.zeros((12,12))
1672    k = 0
1673    for i in range(12):
1674        for j in range(i,12):
1675            vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j]
1676            vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j]
1677            k += 1
1678    M = np.inner(vcov,DYI.T)
1679    WYI = np.sum(M*DYI,axis=0)
1680    WYI = np.where(WYI>0.,WYI,0.)
1681    return YI,WYI
1682   
1683#testing data
1684NeedTestData = True
1685def TestData():
1686    'needs a doc string'
1687#    global NeedTestData
1688    NeedTestData = False
1689    global bakType
1690    bakType = 'chebyschev'
1691    global xdata
1692    xdata = np.linspace(4.0,40.0,36000)
1693    global parmDict0
1694    parmDict0 = {
1695        'pos0':5.6964,'int0':8835.8,'sig0':1.0,'gam0':1.0,
1696        'pos1':11.4074,'int1':3922.3,'sig1':1.0,'gam1':1.0,
1697        'pos2':20.6426,'int2':1573.7,'sig2':1.0,'gam2':1.0,
1698        'pos3':26.9568,'int3':925.1,'sig3':1.0,'gam3':1.0,
1699        'U':1.163,'V':-0.605,'W':0.093,'X':0.0,'Y':2.183,'SH/L':0.002,
1700        'Back0':5.384,'Back1':-0.015,'Back2':.004,
1701        }
1702    global parmDict1
1703    parmDict1 = {
1704        'pos0':13.4924,'int0':48697.6,'sig0':1.0,'gam0':1.0,
1705        'pos1':23.4360,'int1':43685.5,'sig1':1.0,'gam1':1.0,
1706        'pos2':27.1152,'int2':123712.6,'sig2':1.0,'gam2':1.0,
1707        'pos3':33.7196,'int3':65349.4,'sig3':1.0,'gam3':1.0,
1708        'pos4':36.1119,'int4':115829.8,'sig4':1.0,'gam4':1.0,
1709        'pos5':39.0122,'int5':6916.9,'sig5':1.0,'gam5':1.0,
1710        'U':22.75,'V':-17.596,'W':10.594,'X':1.577,'Y':5.778,'SH/L':0.002,
1711        'Back0':36.897,'Back1':-0.508,'Back2':.006,
1712        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1713        }
1714    global parmDict2
1715    parmDict2 = {
1716        'pos0':5.7,'int0':1000.0,'sig0':0.5,'gam0':0.5,
1717        'U':2.,'V':-2.,'W':5.,'X':0.5,'Y':0.5,'SH/L':0.02,
1718        'Back0':5.,'Back1':-0.02,'Back2':.004,
1719#        'Lam1':1.540500,'Lam2':1.544300,'I(L2)/I(L1)':0.5,
1720        }
1721    global varyList
1722    varyList = []
1723
1724def test0():
1725    if NeedTestData: TestData()
1726    msg = 'test '
1727    gplot = plotter.add('FCJ-Voigt, 11BM').gca()
1728    gplot.plot(xdata,getBackground('',parmDict0,bakType,'PXC',xdata)[0])   
1729    gplot.plot(xdata,getPeakProfile(parmDict0,xdata,varyList,bakType))
1730    fplot = plotter.add('FCJ-Voigt, Ka1+2').gca()
1731    fplot.plot(xdata,getBackground('',parmDict1,bakType,'PXC',xdata)[0])   
1732    fplot.plot(xdata,getPeakProfile(parmDict1,xdata,varyList,bakType))
1733   
1734def test1():
1735    if NeedTestData: TestData()
1736    time0 = time.time()
1737    for i in range(100):
1738        y = getPeakProfile(parmDict1,xdata,varyList,bakType)
1739    print '100+6*Ka1-2 peaks=1200 peaks',time.time()-time0
1740   
1741def test2(name,delt):
1742    if NeedTestData: TestData()
1743    varyList = [name,]
1744    xdata = np.linspace(5.6,5.8,400)
1745    hplot = plotter.add('derivatives test for '+name).gca()
1746    hplot.plot(xdata,getPeakProfileDerv(parmDict2,xdata,varyList,bakType)[0])
1747    y0 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1748    parmDict2[name] += delt
1749    y1 = getPeakProfile(parmDict2,xdata,varyList,bakType)
1750    hplot.plot(xdata,(y1-y0)/delt,'r+')
1751   
1752def test3(name,delt):
1753    if NeedTestData: TestData()
1754    names = ['pos','sig','gam','shl']
1755    idx = names.index(name)
1756    myDict = {'pos':parmDict2['pos0'],'sig':parmDict2['sig0'],'gam':parmDict2['gam0'],'shl':parmDict2['SH/L']}
1757    xdata = np.linspace(5.6,5.8,800)
1758    dx = xdata[1]-xdata[0]
1759    hplot = plotter.add('derivatives test for '+name).gca()
1760    hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1])
1761    y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1762    myDict[name] += delt
1763    y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)
1764    hplot.plot(xdata,(y1-y0)/delt,'r+')
1765
1766if __name__ == '__main__':
1767    import GSASIItestplot as plot
1768    global plotter
1769    plotter = plot.PlotNotebook()
1770#    test0()
1771#    for name in ['int0','pos0','sig0','gam0','U','V','W','X','Y','SH/L','I(L2)/I(L1)']:
1772    for name,shft in [['int0',0.1],['pos0',0.0001],['sig0',0.01],['gam0',0.00001],
1773        ['U',0.1],['V',0.01],['W',0.01],['X',0.0001],['Y',0.0001],['SH/L',0.00005]]:
1774        test2(name,shft)
1775    for name,shft in [['pos',0.0001],['sig',0.01],['gam',0.0001],['shl',0.00005]]:
1776        test3(name,shft)
1777    print "OK"
1778    plotter.StartEventLoop()
Note: See TracBrowser for help on using the repository browser.