source: trunk/GSASIIpwd.py @ 1572

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

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