source: trunk/GSASIIpwd.py @ 1571

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

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