source: trunk/GSASIIpwd.py @ 1637

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

fixes to peak indexing routines; new Inst argument for getHKLpeak messed up indexing routines. Now optional at end of arguments.
some more work on SS constraints.

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